// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
// avtBlueprintFileFormat.C
// ************************************************************************* //

#include "avtBlueprintFileFormat.h"

//-----------------------------------------------------------------------------
// visit includes
//-----------------------------------------------------------------------------
#include "avtDatabaseMetaData.h"
#include "avtResolutionSelection.h"

#include "avtMaterial.h"
#include "avtSpecies.h"
#include "avtMixedVariable.h"
#include "avtVariableCache.h"
#include "DBOptionsAttributes.h"

#include "DebugStream.h"
#include "StringHelpers.h"
#include "TimingsManager.h"

#include "Expression.h"
#include "FileFunctions.h"
#include "InvalidVariableException.h"
#include "InvalidFilesException.h"
#include "UnexpectedValueException.h"

//-----------------------------------------------------------------------------
// vtk includes
//-----------------------------------------------------------------------------
#include <vtkUnstructuredGrid.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItUtility.h>

//-----------------------------------------------------------------------------
// std lib includes
//-----------------------------------------------------------------------------
#include <string>

//-----------------------------------------------------------------------------
// conduit includes
//-----------------------------------------------------------------------------
#include "conduit.hpp"
#include "conduit_relay.hpp"
#include "conduit_relay_io_hdf5.hpp"
#include "conduit_blueprint.hpp"
#include "conduit_blueprint_mesh_utils.hpp"
#include "conduit_fmt/conduit_fmt.h"

//-----------------------------------------------------------------------------
// mfem includes
//-----------------------------------------------------------------------------
#include "mfem.hpp"


#ifdef PARALLEL
#include <mpi.h>
#include <avtParallel.h>
#include "conduit_relay_mpi.hpp"
#endif


//-----------------------------------------------------------------------------
// bp visit plugin includes
//-----------------------------------------------------------------------------
#include "avtBlueprintLogging.h"
#include "avtBlueprintTreeCache.h"

#include <vtkDataArray.h>
#include <algorithm>

using std::string;
using namespace conduit;
using namespace mfem;

const std::string avtBlueprintFileFormat::DISPLAY_NAME("display_name");

// ****************************************************************************
//  Method: avtBlueprintFileFormat::FetchMeshAndTopoNames
//
//  Purpose: Maps the full mesh name, registered with visit to its
//           blueprint mesh and topo components.
//
//  Programmer: cyrush
//  Creation:   Mon Mar  9 12:24:16 PDT 2020
//
//  Modifications:
//
//
// ****************************************************************************
void
avtBlueprintFileFormat::FetchMeshAndTopoNames(const std::string &name_name_full,
                                              std::string &mesh_name,
                                              std::string &topo_name)
{
    string mesh_base = FileFunctions::Basename(name_name_full);

    if(!m_mesh_and_topo_info.has_child(mesh_base))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "Unknown mesh name " << mesh_name);
    }

    mesh_name = m_mesh_and_topo_info[mesh_base]["mesh"].as_string();
    topo_name = m_mesh_and_topo_info[mesh_base]["topo"].as_string();
}


// ****************************************************************************
std::string
sanitize_var_name(const std::string &varname)
{
    string res = StringHelpers::Replace(varname,":","_colon_");
    res = StringHelpers::Replace(res,"[","_lb_");
    res = StringHelpers::Replace(res,"]","_rb_");
    return res;
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat constructor
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
// 
//  Modifications:
//    Justin Privitera, Mon Apr 11 18:20:19 PDT 2022
//    Added query to check MFEM LOR Setting. If it is set to legacy, the 
//    constructor sets `m_new_refine` to false, and if it is set to new,
//    the opposite. It is set to new by default.
// 
//    Justin Privitera, Wed Aug 24 11:18:25 PDT 2022
//    Removed setting of info, warning, and error handlers.
//
//    Brad Whitlock, Mon May 22 16:55:01 PDT 2023
//    I initialized some new members.
// 
//    Justin Privitera, Fri Sep 27 11:51:59 PDT 2024
//    Added specset info.
//
// ****************************************************************************
avtBlueprintFileFormat::avtBlueprintFileFormat(const char *filename, DBOptionsAttributes *opts)
    : avtSTMDFileFormat(&filename, 1),
      m_root_node(),
      m_protocol(""),
      m_tree_cache(NULL),
      m_selected_lod(0),
      m_mesh_and_topo_info(),
      m_matset_info(),
      m_specset_info(),
      m_mfem_mesh_map(),
      m_mfem_material_map(),
      m_new_refine(true)
{
    if (opts->GetEnum("MFEM LOR Setting") == 0)
    {
        // legacy LOR was requested
        m_new_refine = false;
    }
    else
    {
        // new LOR was requested
        m_new_refine = true;
    }    

    m_tree_cache = new avtBlueprintTreeCache();
}

avtBlueprintFileFormat::~avtBlueprintFileFormat()
{
    delete m_tree_cache;
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
// 
//  Modifications:
//     Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//     Encased in try-catch block.
//
// ****************************************************************************
void
avtBlueprintFileFormat::FreeUpResources(void)
{
    try
    {
        m_tree_cache->Release();
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::FreeUpResources: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}

// ****************************************************************************
//  Method: GenerateOffsetsForPolytopalMesh
//
//  Purpose:
//      Polyhedral meshes require offsets to verify. We can generate them
//      before verification takes place.
//
//  Arguments:
//      data        The mesh domain we wish to add offsets to.
// 
//  Notes:
//      We cannot assume much about the data we are working with, as it has
//      not passed verification yet. If it has the structure we expect, we can
//      assume that the passed data node consists of a single coordset and
//      topology, as ReadBlueprintMesh() is called directly before this 
//      function and should return a node with a single coordset and topology.
//
//  Programmer: Justin Privitera
//  Creation:   10/25/23
//
//  Modifications:
//
// ****************************************************************************
void
GenerateOffsetsForPolytopalMesh(Node &data)
{
    // using short circuit evaluation, so these statements will only execute
    // until we hit the first one that is false
    if (data.has_child("topologies") && 
        data["topologies"].number_of_children() > 0 && 
        data["topologies"][0].has_child("type") && 
        data["topologies"][0]["type"].as_string() == "unstructured" && 
        data["topologies"][0].has_path("elements/shape") && 
        (data["topologies"][0]["elements/shape"].as_string() == "polyhedral" || 
         data["topologies"][0]["elements/shape"].as_string() == "polygonal") && 
        (! data["topologies"][0].has_path("elements/offsets")))
        blueprint::mesh::utils::topology::unstructured::generate_offsets_inline(data["topologies"][0]);
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadBlueprintMesh
//
//  Purpose:
//      Reads a mesh domain into the `out` conduit Node.
//
//
//  Programmer: Cyrus Harrison and Mark Miller
//  Creation:   Fri Aug 12 13:45:34 PDT 2016
//
//  Modifications:
//    Cyrus Harrison, Mon Mar  9 15:45:17 PDT 2020
//    Use explicit map from registered mesh name to bp mesh and topo names.
//
//    Cyrus Harrison, Tue Mar 10 13:15:26 PDT 2020
//    Support empty mesh case.
//
//    Cyrus Harrison, Tue Dec 13 12:17:17 PST 2022
//    Refactor to pass mesh name to tree cache for part map style index
//    support.
// 
//    Justin Privitera, Wed Oct 25 17:29:07 PDT 2023
//    Call GenerateOffsetsForPolytopalMesh() at the end of this function, 
//    before calling verify.
// 
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
//
// ****************************************************************************
void
avtBlueprintFileFormat::ReadBlueprintMesh(int domain,
                                          const string &abs_meshname,
                                          Node &out)
{
    BP_PLUGIN_INFO("ReadBlueprintMesh: " << abs_meshname
                    << " [domain " << domain << "]");
    string mesh_name;
    string topo_name;
    FetchMeshAndTopoNames(std::string(abs_meshname),
                          mesh_name,
                          topo_name);

    BP_PLUGIN_INFO("mesh name: " << mesh_name);
    BP_PLUGIN_INFO("topo name: " << topo_name);

    if (!m_root_node["blueprint_index"].has_child(mesh_name))
    {
        EXCEPTION1(InvalidVariableException, mesh_name);
    }

    const Node &bp_index_mesh_node = m_root_node["blueprint_index"][mesh_name];
    const Node &bp_index_topo_node = bp_index_mesh_node["topologies"][topo_name];

    string coordset_name = bp_index_topo_node["coordset"].as_string();
    string topo_path     = bp_index_topo_node["path"].as_string();
    string coords_path   = bp_index_mesh_node["coordsets"][coordset_name]["path"].as_string();

    BP_PLUGIN_INFO("coordset path " << coords_path);

    // if we can't fetch the coordset, assume we have an empty domain
    // and skip if that is the case
    try
    {
        m_tree_cache->FetchBlueprintTree(domain,
                                         mesh_name,
                                         coords_path,
                                         out["coordsets"][coordset_name]);
    }
    catch(InvalidVariableException const&)
    {
        BP_PLUGIN_WARNING("failed to load conduit coordset for "
                           << abs_meshname << " [domain "<< domain << "]"
                           << " -- skipping mesh for this domain");
        // if something went wrong, reset the output node to
        // signal the read failed, and return.
        out.reset();
        return;
    }

    BP_PLUGIN_INFO("topology path " << topo_path);

    // if we can't fetch the topo, assume we have an empty domain
    // and skip if that is the case
    try
    {
        m_tree_cache->FetchBlueprintTree(domain,
                                         mesh_name,
                                         topo_path,
                                         out["topologies"][topo_name]);
    }
    catch(InvalidVariableException const&)
    {
        BP_PLUGIN_WARNING("failed to load conduit topo for "
                           << abs_meshname << " [domain "<< domain << "]"
                           << " -- skipping mesh for this domain");
        // if something went wrong, reset the output node to
        // signal the read failed, and return.
        out.reset();
        return;
    }


    BP_PLUGIN_INFO("GetMesh: done loading conduit data for "
                    << abs_meshname << " [domain " <<domain << "]");

    // check for mfem case
    Node &topo_data_node = out["topologies"][topo_name];


    bool has_bndry_topo = topo_data_node.has_child("boundary_topology");
    string bndry_topo_name;

    // if we have an mfem mesh, load extra data necessary to construct
    // a mfem mesh object
    if(has_bndry_topo)
    {
        BP_PLUGIN_INFO(abs_meshname << " has a boundary topology");
        // mfem case, we need to fetch the boundary topo
        bndry_topo_name = topo_data_node["boundary_topology"].as_string();

        BP_PLUGIN_INFO("boundary_topology is named: " << bndry_topo_name);
        if(bp_index_mesh_node["topologies"].has_child(bndry_topo_name))
        {
            const Node &bp_index_boundary_topo_node =bp_index_mesh_node["topologies"][bndry_topo_name];

            string bnd_topo_path =bp_index_boundary_topo_node["path"].as_string();

            string bnd_topo_coordset_name = bp_index_boundary_topo_node["coordset"].as_string();


            // sanity check:
            // make sure the boundary topo is defined on the same coordset as the mesh
            if(bnd_topo_coordset_name != coordset_name)
            {
                BP_PLUGIN_WARNING("warning: boundary topology "
                                   << "'" << bndry_topo_name << "'"
                                   << "coordset (" << bnd_topo_coordset_name << ")"
                                   << " does not match main topology coordset (" << coordset_name << ")");

            }


            BP_PLUGIN_INFO("boundary topology path " << bnd_topo_path);
            m_tree_cache->FetchBlueprintTree(domain,
                                             mesh_name,
                                             bnd_topo_path,
                                             out["topologies"][bndry_topo_name]);
        }
        else
        {

            BP_PLUGIN_WARNING("boundary_topology '"
                               << bndry_topo_name
                               <<  "' not found in blueprint index");
        }

    }

    bool has_grid_func  = topo_data_node.has_child("grid_function");

    if(has_grid_func)
    {
        BP_PLUGIN_INFO(abs_meshname << " has a mfem grid fuction");
        // mfem case, we need to fetch the grid func
        string gf_name = topo_data_node["grid_function"].as_string();
        BP_PLUGIN_INFO("field for grid_function for topology is named: " << gf_name);
        if(!bp_index_mesh_node["fields"].has_child(gf_name))
        {
            BP_PLUGIN_WARNING("grid_function '" << gf_name << "' for topology not found in fields");
        }
        else
        {
           string gf_path = bp_index_mesh_node["fields"][gf_name]["path"].as_string();

           BP_PLUGIN_INFO("grid function path " << gf_path);
           m_tree_cache->FetchBlueprintTree(domain,
                                            mesh_name,
                                            gf_path,
                                            out["fields"][gf_name]);
        }
    }

    // to construct an mfem mesh object,
    // we may also need the mesh and boundary attribute fields
    // here we look for them by name,
    // in the future, we could have them explicitly listed

    if( bp_index_mesh_node.has_child("fields") )
    {
        const Node &bp_index_fields_node = bp_index_mesh_node["fields"];
        NodeConstIterator flds_itr = bp_index_fields_node.children();

        string mesh_att_path;
        string bndry_att_path;

        while(flds_itr.has_next() &&
              // keep going until we have found the path for both
              // the mesh and boundary att fields
              (mesh_att_path.empty() || bndry_att_path.empty()) )
        {
            const Node &fld = flds_itr.next();
            string cld_name = flds_itr.name();
            // see if the name contains "attribute"
            std::size_t att_loc = cld_name.find("attribute");

            // check if it is defined on the mesh topo or boundary topo
            if(att_loc != std::string::npos)
            {
                // we found a contender, see if its topo is the main mesh topo or
                // the boundary topo
                string fld_topo_name = fld["topology"].as_string();

                if(fld_topo_name == topo_name)
                {
                    mesh_att_path = bp_index_fields_node[cld_name]["path"].as_string();
                }
                else if(fld_topo_name == bndry_topo_name)
                {
                    bndry_att_path = bp_index_fields_node[cld_name]["path"].as_string();
                }

            }
        }


        if(!mesh_att_path.empty())
        {
            // load the data for the mesh att tree
            m_tree_cache->FetchBlueprintTree(domain,
                                             mesh_name,
                                             mesh_att_path,
                                             out["fields/mesh_attribute"]);
        }

        if(!bndry_att_path.empty())
        {
            // load the data for the bndry att tree
            m_tree_cache->FetchBlueprintTree(domain,
                                             mesh_name,
                                             bndry_att_path,
                                             out["fields/boundary_attribute"]);
        }
    }

    // before we verify, we want to generate offsets for unstructured polytopal
    // meshes if they do not already exist
    GenerateOffsetsForPolytopalMesh(out);

}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetBlueprintIndexForField
//
//  Purpose:
//    Returns a node in the Blueprint index for the field that either has the
//    supplied name or has a DISPLAY_NAME that matches the supplied name.
//
//  Programmer: Brad Whitlock
//  Creation:   Tue Jul 18 15:57:00 PDT 2023
//
//  Modifications:
//
// ****************************************************************************

const Node *
avtBlueprintFileFormat::GetBlueprintIndexForField(const Node &fields,
    const string &abs_varname) const
{
    const char *mName = "GetBlueprintIndexForField";
    BP_PLUGIN_INFO(mName << ": " << abs_varname);
    const Node *retval = nullptr;
    string varname = FileFunctions::Basename(abs_varname);

    if (fields.has_child(varname))
    {
        BP_PLUGIN_INFO(mName << ": found " << varname);
        retval = fields.fetch_ptr(varname);
    }
    else
    {
        // We may have been given a display name for the variable as opposed
        // to its real name. See if we can find that.
        BP_PLUGIN_INFO(mName << ": checking fields for a display name.");
        for(conduit::index_t i = 0; i < fields.number_of_children(); i++)
        {
            const Node *f = fields.child_ptr(i);
            BP_PLUGIN_INFO(mName << ": checking field " << f->name() << ".");
            if(f->has_child(DISPLAY_NAME))
            {
                std::string dn = f->fetch_existing(DISPLAY_NAME).as_string();
                if(dn == abs_varname)
                {
                    BP_PLUGIN_INFO(mName << ": " << f->name() << " matched based on "
                                   << DISPLAY_NAME << " " << dn);
                    retval = f;
                    break;
                }
            }
        }

        // only throw an error if element_coloring is not in the name
        // element_coloring won't be in the index, its automatic.
        if(retval == nullptr &&
           varname.find("element_coloring") == std::string::npos)
        {
            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                 "field " << varname << " not found in blueprint index");
        }
    }
    return retval;
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadBlueprintField
//
//  Purpose:
//      Reads a field for the given domain into the `out` conduit Node.
//
//
//  Programmer: Cyrus Harrison and Mark Miller
//  Creation:   Fri Aug 12 13:45:34 PDT 2016
//
//  Modifications:
//    Cyrus Harrison, Mon Mar  9 15:45:17 PDT 2020
//    Use explicit map from registered mesh name to bp mesh and topo names.
//
//    Cyrus Harrison, Wed Mar 11 10:42:22 PDT 2020
//    Allow empty domains.
// 
//    Justin Privitera, Wed Mar 22 16:09:52 PDT 2023
//    Handle the 1D curve case.
// 
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
//
//    Brad Whitlock, Tue Jul 18 16:02:13 PDT 2023
//    Locate the index node for the variable using GetBlueprintIndexField.
//
// ****************************************************************************

void
avtBlueprintFileFormat::ReadBlueprintField(int domain,
                                           const string &abs_varname,
                                           Node &out)
{
    BP_PLUGIN_INFO("ReadBlueprintField: " << abs_varname << " [domain " << domain << "]");

    string abs_varname_str(abs_varname);
    // replace colons, etc
    abs_varname_str = sanitize_var_name(abs_varname_str);
    string abs_meshname = metadata->MeshForVar(abs_varname_str);

    if (m_curve_names.find(abs_meshname) != m_curve_names.end())
    {
        // This means that abs_meshname is of the form:
        // "mesh_name/field_name", so we need to separate it
        abs_meshname = abs_meshname.substr(0, abs_meshname.find('/'));
    }

    BP_PLUGIN_INFO("field " << abs_varname << " is defined on mesh " << abs_meshname);

    string mesh_name;
    string topo_name;
    FetchMeshAndTopoNames(abs_meshname,
                          mesh_name,
                          topo_name);

    BP_PLUGIN_INFO("mesh name: " << mesh_name);
    BP_PLUGIN_INFO("topo name: " << topo_name);

    if (!m_root_node["blueprint_index"].has_child(mesh_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "mesh " << mesh_name << " not found in blueprint index");
    }

    // Look up the index node for the field.
    const Node &fields = m_root_node["blueprint_index"][mesh_name]["fields"];
    const Node *bp_index_field = GetBlueprintIndexForField(fields, abs_varname);
    std::string data_path;
    if(bp_index_field != nullptr)
    {
        BP_PLUGIN_INFO(bp_index_field->to_yaml());
        data_path    = bp_index_field->fetch_existing("path").as_string();
    }

    try
    {
        m_tree_cache->FetchBlueprintTree(domain,
                                         mesh_name,
                                         data_path,
                                         out);
        BP_PLUGIN_INFO("done loading conduit data for " << abs_varname << " [domain "<< domain << "]" );
    }
    catch(InvalidVariableException const&)
    {
        BP_PLUGIN_WARNING("failed to load conduit data for "
                           << abs_varname << " [domain "<< domain << "]"
                           << " -- skipping field for this domain");
        // if something went wrong, reset the output node to
        // signal the read failed, and return.
        out.reset();
    }
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::DetectHOMaterial
//
//  Purpose:
//    Determines whether there is a HO material.
//
//
//  Programmer: Brad Whitlock
//  Creation:   Mon May 22 16:51:12 PDT 2023
//
//  Modifications:
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
//
//     Cyrus Harrison Thu Nov  2 11:23:40 PDT 2023
//     Added new detection case for "volume_fraction_" prefix.
//
// ****************************************************************************

bool
avtBlueprintFileFormat::DetectHOMaterial(const std::string &mesh_name,
    const std::string &topo_name,
    const std::vector<std::string> &matNames,
    std::map<std::string, std::string> &matFields, std::string &freeMatName) const
{
    bool HOmaterials = false;
    bool axom_vol_frac_convention = false;
    freeMatName.clear();
    if (m_root_node["blueprint_index"][mesh_name].has_child("fields"))
    {
        const conduit::Node &n_fields = m_root_node["blueprint_index"][mesh_name]["fields"];
        // Look for Axom convention.
        std::vector<std::string> prefix_list = {"vol_frac_","volume_fraction_"};

        for(size_t i = 0; i < matNames.size(); i++)
        {
            const auto &matname = matNames[i];
            for(size_t prefix_idx = 0; prefix_idx < prefix_list.size(); prefix_idx++)
            {
                std::string fieldName(prefix_list[prefix_idx] + matname);
                if(n_fields.has_child(fieldName))
                {
                    if(prefix_list[prefix_idx] == "vol_frac_")
                    {
                        axom_vol_frac_convention = true;
                    }

                    const conduit::Node &f = n_fields.fetch_existing(fieldName);
                    if(f.has_child("basis") &&
                    f.has_child("topology") &&
                    f.has_child("number_of_components"))
                    {
                        std::string tname = f["topology"].as_string();
                        int nc = f["number_of_components"].to_int();
                        if(tname == topo_name && nc == 1)
                        {
                            matFields[matname] = fieldName;
                        }
                    }
                }
            }
        }
        // If all of the material names had a matching HO field for the volume
        // fractions then our material is made up of HO fields.
        HOmaterials = matFields.size() == matNames.size();

        // See whether a free material needs to be created. Use Axom convention.
        if(axom_vol_frac_convention)
        {
            // See whether a free material needs to be created. Use Axom convention.
            const std::string free_mat_name("vol_frac_free");
            bool make_free_mat = !m_root_node["blueprint_index"][mesh_name]["fields"].has_child(free_mat_name);
            if(make_free_mat)
                freeMatName = "free";
        }
    }
    return HOmaterials;
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadBlueprintMatset
//
//  Purpose:
//      Reads matset info for the given domain into the `out` conduit Node.
//
//
//  Programmer: Cyrus Harrison
//  Creation:   Wed Dec  9 13:02:46 PST 2020
//
//  Modifications:
//    Brad Whitlock, Mon May 22 16:51:12 PDT 2023
//    I added code to treat HO materials specially since we want them to be
//    refined according to the selected level of detail (m_selected_lod).
// 
//    Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//    Fixed warnings.
//
//    Brad Whitlock, Wed Jul 19 13:56:42 PDT 2023
//    I added display_name support.
//
// ****************************************************************************

void
avtBlueprintFileFormat::ReadBlueprintMatset(int domain,
                                            const string &abs_matsetname,
                                            Node &out)
{
    BP_PLUGIN_INFO("ReadBlueprintMatsetVolFracs: " << abs_matsetname
                   << " [domain " << domain << "]");

    // replace colons, etc
    const std::string abs_matsetname_str = sanitize_var_name(abs_matsetname);

    // we need to know what mesh name this matset is associated with
    const Node &mset_info = m_matset_info[abs_matsetname_str];
    
    const string abs_meshname = mset_info["full_mesh_name"].as_string();

    BP_PLUGIN_INFO("matset " << abs_matsetname << " is defined on mesh " << abs_meshname);

    string mesh_name = mset_info["mesh_name"].as_string();
    string topo_name = mset_info["topo_name"].as_string();
    string matset_name = mset_info["matset_name"].as_string();
    const Node &n_mat_names = mset_info["matnames"];

    BP_PLUGIN_INFO("mesh name: " << mesh_name);
    BP_PLUGIN_INFO("topo name: " << topo_name);
    BP_PLUGIN_INFO("matset name: " << matset_name);
    BP_PLUGIN_INFO("matnames: " << n_mat_names.to_yaml());

    if (!m_root_node["blueprint_index"].has_child(mesh_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "mesh " << mesh_name << " not found in blueprint index");
    }

    if (!m_root_node["blueprint_index"][mesh_name]["matsets"].has_child(matset_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "matset " << matset_name << " not found in blueprint index");
    }

    const Node &bp_index_matset = m_root_node["blueprint_index"][mesh_name]["matsets"][matset_name];
    BP_PLUGIN_INFO(bp_index_matset.to_yaml());

    const string data_path = bp_index_matset["path"].as_string();

    // See whether the materials in the index correspond to HO fields.
    std::vector<std::string> matNames;
    for(conduit::index_t i = 0; i < n_mat_names.number_of_children(); i++)
        matNames.push_back(n_mat_names[i].name());
    std::map<std::string, std::string> matFields;
    std::string freeMatName;
    if(DetectHOMaterial(mesh_name, topo_name, matNames, matFields, freeMatName))
    {
        const std::string line("=============================================================================");
        BP_PLUGIN_INFO(line << endl
                            << " Start Reading HO material " << abs_matsetname_str << endl
                            << line);

        // See whether a free material needs to be created.
        bool make_free_mat = !freeMatName.empty();
        std::vector<float> freevf;

        conduit::Node &vf = out["volume_fractions"];
        conduit::Node &mn = out["matnames"];
        int idx = 0;
        for(auto it = matFields.begin(); it != matFields.end(); it++)
        {
            // Make a variable name for the volume fraction field and read it.
            std::string matVar = mesh_name + "_" + topo_name + "/" + it->second;
            // If there is a display_name, use that to request the data.
            const Node &bp_idx_fields = m_root_node["blueprint_index"][mesh_name]["fields"];
            const Node *bpi = GetBlueprintIndexForField(bp_idx_fields, it->second);
            if(bpi != nullptr && bpi->has_child(DISPLAY_NAME))
                matVar = bpi->fetch_existing(DISPLAY_NAME).as_string();
            vtkDataArray *da = GetVar(domain, matVar.c_str());

            // Save the material field as float into the new out node.
            conduit::Node &f = vf[it->first];
            auto nzones = static_cast<vtkIdType>(da->GetNumberOfTuples());
            f.set(conduit::DataType::float32(nzones));
            float *fptr = f.as_float32_ptr();
            for(vtkIdType zi = 0; zi < nzones; ++zi)
                fptr[zi] = static_cast<float>(da->GetTuple1(zi));
            da->Delete();

            // If we need to make a free material, do it.
            if(make_free_mat)
            {
                if(freevf.empty())
                    freevf.resize(static_cast<size_t>(nzones), 1.f);
                // Now, subtract the current vf from the free mat.
                for(vtkIdType zi = 0; zi < nzones; ++zi)
                    freevf[static_cast<size_t>(zi)] -= fptr[zi];
            }

            // Add the material name to the list.
            mn[it->first] = idx++;
        }
        // Append the free material.
        if(make_free_mat)
        {
            // See whether any zones have sufficient free material.
            auto it = std::find_if(freevf.begin(), freevf.end(), [](float value)
            {
                constexpr float SUFFICIENT_MATERIAL = 1.e-6f;
                return (1.f - value) > SUFFICIENT_MATERIAL;
            });
            if(it != freevf.end())
            {
                vf[freeMatName].set(freevf);
                mn[freeMatName] = idx++;
            }
        }
        out["topology"] = mesh_name;

        // Save the LOD used to make the material. This way we know if we need
        // to purge it the next time we ask for the mesh.
        BP_PLUGIN_INFO("Record that material " << abs_matsetname
                       << " is at LOD " << (m_selected_lod + 1));
        m_mfem_material_map[topo_name] = std::make_pair(abs_matsetname_str, m_selected_lod + 1);

        BP_PLUGIN_INFO(line << endl
                            << " Done Reading HO material " << abs_matsetname_str << endl
                            << line);
    }
    else
    {
        // Non-HO path.
        try
        {
            m_tree_cache->FetchBlueprintTree(domain,
                                             mesh_name,
                                             data_path,
                                             out);

            BP_PLUGIN_INFO("done loading conduit data for " 
                            << abs_matsetname << " [domain "<< domain << "]" );
        }
        catch(InvalidVariableException const&)
        {
            BP_PLUGIN_WARNING("failed to load conduit data for "
                               << abs_matsetname << " [domain "<< domain << "]"
                               << " -- skipping field for this domain");
            // if something went wrong, reset the output node to
            // signal the read failed, and return.
            out.reset();
        }

        // provide material_map
        out["matnames"] = n_mat_names;
    }
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadBlueprintSpecset
//
//  Purpose:
//      Reads specset info for the given domain into the `out` conduit Node.
//
//
//  Programmer: Justin Privitera
//  Creation:   09/27/24
//
//  Modifications:
//
// ****************************************************************************

void
avtBlueprintFileFormat::ReadBlueprintSpecset(int domain,
                                            const string &specset_name,
                                            Node &out)
{
    BP_PLUGIN_INFO("ReadBlueprintSpecsetMassFracs: " << specset_name
                   << " [domain " << domain << "]");

    // replace colons, etc
    const std::string specset_name_str = sanitize_var_name(specset_name);

    // we need to know what mesh name this specset is associated with
    const Node &specset_info = m_specset_info[specset_name_str];
    const std::string abs_meshname = specset_info["full_mesh_name"].as_string();
    const std::string matset_name = specset_info["matset_name"].as_string();
    const std::string bp_specset_name = specset_info["specset_name"].as_string();

    BP_PLUGIN_INFO("specset " << specset_name << " is defined on mesh " << abs_meshname);

    const std::string mesh_name = specset_info["mesh_name"].as_string();
    const std::string topo_name = specset_info["topo_name"].as_string();
    const Node &n_spec_names = specset_info["species"];

    BP_PLUGIN_INFO("mesh name: " << mesh_name);
    BP_PLUGIN_INFO("topo name: " << topo_name);
    BP_PLUGIN_INFO("matset name: " << matset_name);
    BP_PLUGIN_INFO("specset name: " << bp_specset_name);
    BP_PLUGIN_INFO("specnames: " << n_spec_names.to_yaml());

    if (!m_root_node["blueprint_index"].has_child(mesh_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "mesh " << mesh_name << " not found in blueprint index");
    }

    if (!m_root_node["blueprint_index"][mesh_name]["matsets"].has_child(matset_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "matset " << matset_name << " not found in blueprint index");
    }

    if (!m_root_node["blueprint_index"][mesh_name]["specsets"].has_child(bp_specset_name))
    {
        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                             "specset " << bp_specset_name << " not found in blueprint index");
    }

    const Node &bp_index_specset = m_root_node["blueprint_index"][mesh_name]["specsets"][bp_specset_name];
    BP_PLUGIN_INFO(bp_index_specset.to_yaml());

    const string data_path = bp_index_specset["path"].as_string();

    try
    {
        m_tree_cache->FetchBlueprintTree(domain,
                                         mesh_name,
                                         data_path,
                                         out);

        BP_PLUGIN_INFO("done loading conduit data for " 
                        << specset_name_str << " [domain "<< domain << "]" );
    }
    catch(InvalidVariableException const&)
    {
        BP_PLUGIN_WARNING("failed to load conduit data for "
                           << specset_name_str << " [domain "<< domain << "]"
                           << " -- skipping field for this domain");
        // if something went wrong, reset the output node to
        // signal the read failed, and return.
        out.reset();
    }
}


// ****************************************************************************
// helper method used to add the meta data for a blueprint mesh.
// ****************************************************************************
//  Modifications:
//    Cyrus Harrison, Mon Mar  9 15:45:17 PDT 2020
//    Use explicit map from registered mesh name to bp mesh and topo names.
//
//    Justin Privitera, Mon 07 Mar 2022 03:08:01 PM PST
//    Added case for implicit points topology.
// 
//    Justin Privitera, Fri Jul 22 16:10:43 PDT 2022
//    Added check for mfem zone centered field.
// 
//    Justin Privitera, Wed Oct 19 15:03:26 PDT 2022
//    Added logic to determine nodal vs zonal association for mfem grid 
//    functions.
// 
//    Justin Privitera, Wed Mar 22 16:09:52 PDT 2023
//    Handle 1D curve case.
// 
//    Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//    Fixed warnings.
// 
//    Justin Privitera, Wed Oct 30 14:18:31 PDT 2024
//    Determine if MFEM meshes could be periodic (all periodic are L2 but not 
//    all L2 are periodic) and mark fields for those meshes as nodal because
//    they will use legacy LOR down the line which makes all fields nodal.
//
// ****************************************************************************
void
avtBlueprintFileFormat::AddBlueprintMeshAndFieldMetadata(avtDatabaseMetaData *md,
                                                         string const &mesh_name,
                                                         const Node &n_mesh_info)
{


    Node verify_info;
    if(!blueprint::mesh::index::verify(n_mesh_info,verify_info))
    {
        BP_PLUGIN_INFO("Skipping mesh named \"" << mesh_name << "\"" << endl
                       << "blueprint::mesh::index::verify failed " << endl
                       << verify_info.to_yaml());
        return;
    }

    BP_PLUGIN_INFO("Adding mesh named \"" << mesh_name << "\"");

    const Node &n_topos = n_mesh_info["topologies"];

    // holds the names of the topos we found
    // so we can do a sanity check when we add
    // fields
    std::map<std::string,int> topo_dims;

    // a map from topology names to grid function names
    // ... if any exist
    std::map<std::string, std::string> topo_names_to_gf_names;

    //
    // loop over topologies
    //
    NodeConstIterator topos_itr = n_topos.children();

    while(topos_itr.has_next())
    {
        avtMeshType mt = AVT_UNKNOWN_MESH;

        const Node &n_topo = topos_itr.next();
        string topo_name = topos_itr.name();

        // add info that maps the full mesh name registered with VisIt
        // to the blueprint the mesh and topology names
        string mesh_topo_name = mesh_name + "_" + topo_name;

        m_mesh_and_topo_info[mesh_topo_name]["mesh"] = mesh_name;
        m_mesh_and_topo_info[mesh_topo_name]["topo"] = topo_name;

        bool is_mfem_mesh = false;

        if(n_topo.has_child("grid_function"))
        {
            BP_PLUGIN_INFO(mesh_topo_name << " is an mfem mesh");
            is_mfem_mesh = true;
            topo_names_to_gf_names[topo_name] = n_topo["grid_function"].as_string();
        }

        string coordset_name = n_topo["coordset"].as_string();

        BP_PLUGIN_INFO("topology " << topo_name << " references 'coordset' "
               << coordset_name);

        if (n_topo["type"].as_string() == "uniform")
        {
            mt = AVT_RECTILINEAR_MESH;
            BP_PLUGIN_INFO(mesh_topo_name << " topology is uniform ");

        }
        else if (n_topo["type"].as_string() == "rectilinear")
        {
            mt = AVT_RECTILINEAR_MESH;
            BP_PLUGIN_INFO(mesh_topo_name << " topology is rectilinear ");
        }
        else if (n_topo["type"].as_string() == "structured")
        {
            mt = AVT_CURVILINEAR_MESH;
            BP_PLUGIN_INFO(mesh_topo_name << " topology is structured ");
        }
        else if (n_topo["type"].as_string() == "points")
        {
            mt = AVT_UNSTRUCTURED_MESH;
            BP_PLUGIN_INFO(mesh_topo_name << " topology is unstructured (points) ");
        }
        else if (n_topo["type"].as_string() == "unstructured")
        {
            mt = AVT_UNSTRUCTURED_MESH;
            BP_PLUGIN_INFO(mesh_topo_name << " topology is unstructured ");
        }
        else
        {

            BP_PLUGIN_INFO( "Encountered unknown topology type, \""
                            << n_topo["type"].as_string() << "\"" << endl
                            << "Skipping this mesh for now");
            return;
        }

        //
        // Get number of blocks
        //
        const Node &n_state = n_mesh_info["state"];
        int nblocks = n_state["number_of_domains"].to_int();

        BP_PLUGIN_INFO("number_of_domains: " << nblocks);

        //
        // Get the mesh spatial dimensions
        //

        BP_PLUGIN_INFO("finding coordinate system");

        const Node &n_coordsets = n_mesh_info["coordsets"];
        const Node &n_coords = n_coordsets[coordset_name];

        std::string coord_sys_type = n_coords["coord_system/type"].as_string();

        int ndims = static_cast<int>(n_coords["coord_system/axes"].number_of_children());
        topo_dims[topo_name] = ndims;

        BP_PLUGIN_INFO("coordinate system: "
                        << n_coords["coord_system"].to_yaml()
                        << " (type=" << coord_sys_type
                        << " ndims=" << ndims << ")");

        avtMeshMetaData *mmd = new avtMeshMetaData(mesh_topo_name,
                                                   nblocks,
                                                   0, 0, 0,
                                                   ndims, ndims, mt);

        if(is_mfem_mesh)
        {
            mmd->containsOriginalCells = true;
        }

        if(coord_sys_type == "cylindrical")
        {
            mmd->meshCoordType = AVT_RZ;
            mmd->xLabel = "Z-Axis";
            mmd->yLabel = "R-Axis";
        }
        else
        {
            mmd->meshCoordType = AVT_XY;
        }

        mmd->LODs = 20;
        md->Add(mmd);

        if(is_mfem_mesh)
        {
            // if we have a mfem mesh, add extra element_color variable
            md->Add(new avtScalarMetaData(mesh_topo_name + "/element_coloring",
                                          mesh_topo_name,
                                          AVT_ZONECENT));

            m_mfem_mesh_map[mesh_topo_name] = true;
        }
        else
        {
            m_mfem_mesh_map[mesh_topo_name] = false;
        }
    }

    BP_PLUGIN_INFO("adding field vars for " <<  mesh_name);

    //
    // Now, handle any fields defined for this mesh
    //

    if(n_mesh_info.has_child("fields"))
    {
        // This is a set of topology names that MAY BE periodic.
        // There is no way to check directly, so we are forced to
        // filter all L2 meshes.
        std::set<std::string> periodic_topos;

        // examine mfem basis functions to filter periodic meshes
        for (const auto & topo_gf : topo_names_to_gf_names)
        {
            const std::string topo_name = topo_gf.first;
            const std::string gf_name = topo_gf.second;
            if (n_mesh_info["fields"].has_child(gf_name))
            {
                const Node &n_field = n_mesh_info["fields"][gf_name];
                if (n_field.has_child("basis"))
                {
                    const std::string basis = n_field["basis"].as_string();
                    if (basis.find("L2_") != std::string::npos)
                    {
                        // COULD BE periodic
                        periodic_topos.insert(topo_name);
                    }
                }
            }
        }

        NodeConstIterator fields_itr = n_mesh_info["fields"].children();

        while (fields_itr.has_next())
        {
            const Node &n_field = fields_itr.next();
            string varname = fields_itr.name();
            string var_topo_name = n_field["topology"].as_string();
            string var_mesh_name = mesh_name + "_" + var_topo_name;
            string varname_wmesh = var_mesh_name + "/" + varname;

            // Make the variable name that the user sees.
            string varname_display(varname_wmesh);
            if(n_field.has_child(DISPLAY_NAME))
            {
                varname_display = n_field[DISPLAY_NAME].as_string();
                BP_PLUGIN_INFO("Field \"" << varname
                               << "\" is being exposed as \""
                               << varname_display << "\".");
            }

            if (topo_dims[var_topo_name] == 0)
            {
                BP_PLUGIN_WARNING("Field \"" << varname_display
                                  << "\" defined on unknown topology=\""
                                  << n_field["topology"].as_string());
                continue;
            }

            int ncomps = n_field["number_of_components"].to_int();
            int ndims = topo_dims[var_topo_name];
            
            // 
            // handle centering
            // 
            avtCentering cent = AVT_NODECENT; // default
            if (n_field.has_child("association") &&
                n_field["association"].as_string() == "element")
            {
                cent = AVT_ZONECENT;
            }
            else if (n_field.has_child("basis"))
            {
                // if any of the fields are mfem grid funcs, we may have to
                // treat the mesh as an mfem mesh, even if it lacks a basis func

                m_mfem_mesh_map[var_topo_name] = true;

                if (periodic_topos.count(var_topo_name) > 0)
                {
                    // if this field belongs to a topology that might be a periodic 
                    // mfem mesh then we are always nodal because we are going to
                    // fall back to legacy LOR.
                    cent = AVT_NODECENT;
                }
                else if (m_new_refine) // if new LOR is turned on
                {
                    const std::string basis = n_field["basis"].as_string();
                    // H1 is nodal
                    // L2 is zonal
                    const bool l2 = basis.find("L2_") != std::string::npos;
                    const bool h1 = basis.find("H1_") != std::string::npos;
                    bool node_centered;
                    if (h1 && l2)
                    {
                        BP_PLUGIN_EXCEPTION1(InvalidVariableException, 
                            "AddBlueprintMeshAndFieldMetadata: grid function cannot be both H1 and L2");
                    }
                    else if (!h1 && !l2) // guess
                    {
                        BP_PLUGIN_INFO("WARNING: AddBlueprintMeshAndFieldMetadata: Grid Function is "
                                      "neither H1 nor L2. Guessing nodal association.");
                        node_centered = true; 
                    }
                    else
                    {
                        node_centered = h1 && !l2;
                    }

                    if (!node_centered)
                    {
                        cent = AVT_ZONECENT;
                    }
                }
            }

            // special 1D case
            if (ndims == 1)
            {
                m_curve_names.insert(varname_display);
                avtCurveMetaData *curve = new avtCurveMetaData;
                curve->name = varname_display;
                md->Add(curve);
            }
            else if (ncomps == 1)
                md->Add(new avtScalarMetaData(varname_display, var_mesh_name, cent));
            else if (ndims == 2 && ncomps == 2)
                md->Add(new avtVectorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else if (ndims == 2 && ncomps == 3)
                md->Add(new avtSymmetricTensorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else if (ndims == 2 && ncomps == 4)
                md->Add(new avtTensorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else if (ndims == 3 && ncomps == 3)
                md->Add(new avtVectorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else if (ndims == 3 && ncomps == 6)
                md->Add(new avtSymmetricTensorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else if (ndims == 3 && ncomps == 9)
                md->Add(new avtTensorMetaData(varname_display, var_mesh_name, cent, ncomps));
            else
                md->Add(new avtArrayMetaData(varname_display, var_mesh_name, cent, ncomps));
        }
    }
}


// ****************************************************************************
// helper method used to add materials meta data for a blueprint mesh.
//
// Cyrus Harrison, Tue Dec  8 10:29:21 PST 2020
//
// Modifications:
//   Brad Whitlock, Tue May 23 16:04:30 PDT 2023
//   Added some special handling if the materials are HO.
// 
//   Justin Privitera, Tue Sep 19 11:36:45 PDT 2023
//   No longer assume material ids will be in the range [0, N) and no longer
//   sort the names before sending to avtMaterialMetaData.
// 
//   Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//   Fixed warnings.
// 
//   Justin Privitera, Wed Feb 14 11:37:06 PST 2024
//   Present material ids alongside material names.
// 
//   Justin Privitera, Fri Mar 15 15:56:13 PDT 2024
//   Revert previous change.
//
// ****************************************************************************
void
avtBlueprintFileFormat::AddBlueprintMaterialsMetadata(avtDatabaseMetaData *md,
                                                      string const &mesh_name,
                                                      const Node &n_mesh_info)
{
    if (!n_mesh_info.has_child("matsets"))
    {
        BP_PLUGIN_INFO("Input data file has no matsets.");
        return;
    }
    BP_PLUGIN_INFO("adding materials for " <<  mesh_name);

    NodeConstIterator msets_itr = n_mesh_info["matsets"].children();

    while (msets_itr.has_next())
    {
        const Node &n_mset = msets_itr.next();
        string mset_name = msets_itr.name();

        // the material names are in the "materials" or "material_map" index
        // entries
        if( !n_mset.has_child("materials") && 
            !n_mset.has_child("material_map") )
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " matset index: "
                           << mset_name
                           << " missing `material_map` or `materials`,"
                           << " skipping matset" );
            return;
        }

        // we also need the associated topo
        if (!n_mset.has_child("topology"))
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " matset index: "
                           << mset_name
                           << " missing `topology`,"
                           << " skipping matset" );
            return;
        }

        const std::string topo_name = n_mset["topology"].as_string();
        const string mesh_topo_name = mesh_name + "_" + topo_name;

        const string mesh_matset_name = mesh_topo_name + "_" + mset_name;

        BP_PLUGIN_INFO("adding material set "
                        <<  mesh_topo_name << " " <<  mesh_matset_name);

        if ( n_mset.has_child("material_map") )
        {
            BP_PLUGIN_INFO("material map " << n_mset["material_map"].to_yaml());
        
            NodeConstIterator itr = n_mset["material_map"].children();
            while (itr.has_next())
            {
                const Node &curr_mat = itr.next();
                const int32 mat_id = curr_mat.to_int32();
                const std::string matname = itr.name();
                m_matset_info[mesh_matset_name]["matnames"][matname] = mat_id;
            }
        }
        else // "materials" case, old path
        {
            BP_PLUGIN_INFO("material names " << n_mset["materials"].to_yaml());

            NodeConstIterator itr = n_mset["materials"].children();
            while (itr.has_next())
            {
                itr.next();
                int32 mat_id = static_cast<int32>(itr.index());
                std::string mat_name = itr.name();
                // cache mat names and idx (implied order)
                m_matset_info[mesh_matset_name]["matnames"][mat_name] = mat_id;
            }
        }

        // get matnames vec. No need to sort
        std::vector<string> matnames = m_matset_info[mesh_matset_name]["matnames"].child_names();

        // If the materials were HO then we may need to add a "free" material
        // to the list.
        std::map<std::string, std::string> matFields;
        std::string freeMatName;
        if(DetectHOMaterial(mesh_name, topo_name, matnames, matFields, freeMatName))
        {
            if(!freeMatName.empty())
                matnames.push_back(freeMatName);
        }

        m_matset_info[mesh_matset_name]["full_mesh_name"] = mesh_topo_name;
        m_matset_info[mesh_matset_name]["mesh_name"] = mesh_name;
        m_matset_info[mesh_matset_name]["topo_name"] = topo_name;
        m_matset_info[mesh_matset_name]["matset_name"] = mset_name;
        
        BP_PLUGIN_INFO("Matset Info for "
                       << mesh_matset_name
                       << " : " << m_matset_info[mesh_matset_name].to_yaml())

        avtMaterialMetaData *mmd = new avtMaterialMetaData(mesh_matset_name,
                                                           mesh_topo_name,
                                                           static_cast<int>(matnames.size()),
                                                           matnames);

        mmd->validVariable = true;
        mmd->hideFromGUI = false;
        md->Add(mmd);
    }
}


// ****************************************************************************
// helper method used to add species meta data for a blueprint mesh.
//
// Justin Privitera, Wed Mar 13 16:18:30 PDT 2024
//
// Modifications:
//
// ****************************************************************************
void
avtBlueprintFileFormat::AddBlueprintSpeciesMetadata(avtDatabaseMetaData *md,
                                                    string const &mesh_name,
                                                    const Node &n_mesh_info)
{
    if (!n_mesh_info.has_child("specsets"))
    {
        BP_PLUGIN_INFO("Input data file has no specsets.");
        return;
    }
    BP_PLUGIN_INFO("adding species for " <<  mesh_name);

    NodeConstIterator specsets_itr = n_mesh_info["specsets"].children();

    while (specsets_itr.has_next())
    {
        const Node &n_specset = specsets_itr.next();
        const string specset_name = specsets_itr.name();

        if (!n_specset.has_child("matset"))
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " specset index: "
                           << specset_name
                           << " missing `matset`,"
                           << " skipping specset" );
            return;
        }

        if (!n_specset.has_child("species"))
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " specset index: "
                           << specset_name
                           << " missing `species`,"
                           << " skipping specset" );
            return;
        }

        const std::string matset_name = n_specset["matset"].as_string();

        if (! n_mesh_info.has_path("matsets/" + matset_name))
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " specset index: "
                           << specset_name
                           << " associated matset missing,"
                           << " skipping specset" );
            return;
        }

        if (! n_mesh_info.has_path("matsets/" + matset_name + "/topology"))
        {
            BP_PLUGIN_INFO("mesh: "
                           << mesh_name
                           << " specset index: "
                           << specset_name
                           << " associated matset missing `topology`,"
                           << " skipping specset" );
            return;
        }

        const std::string topo_name = n_mesh_info["matsets"][matset_name]["topology"].as_string();
        const string mesh_topo_name = mesh_name + "_" + topo_name;
        const string mesh_matset_name = mesh_topo_name + "_" + matset_name;
        const string mesh_specset_name = mesh_topo_name + "_" + matset_name + "_" + specset_name;

        BP_PLUGIN_INFO("adding species set "
                        <<  mesh_topo_name << " " <<  mesh_specset_name);

        const int nmat = n_specset["species"].number_of_children();

        std::vector<int> num_species;
        std::vector<std::vector<std::string>> species_names;

        auto mat_itr = n_specset["species"].children();
        while (mat_itr.has_next())
        {
            const Node &mat_val = mat_itr.next();
            num_species.push_back(mat_val.number_of_children());

            std::vector<std::string> spec_names_for_mat;
            auto specie_itr = mat_val.children();
            while (specie_itr.has_next())
            {
                const std::string specname = specie_itr.next().name();
                spec_names_for_mat.push_back(specname);
            }
            species_names.push_back(spec_names_for_mat);
        }

        m_specset_info[mesh_specset_name]["full_mesh_name"] = mesh_topo_name;
        m_specset_info[mesh_specset_name]["mesh_name"] = mesh_name;
        m_specset_info[mesh_specset_name]["topo_name"] = topo_name;
        m_specset_info[mesh_specset_name]["full_material_name"] = mesh_matset_name;
        m_specset_info[mesh_specset_name]["matset_name"] = matset_name;
        m_specset_info[mesh_specset_name]["specset_name"] = specset_name;

        m_specset_info[mesh_specset_name]["nmat"] = nmat;
        m_specset_info[mesh_specset_name]["nmatspec"].set(num_species);
        m_specset_info[mesh_specset_name]["species"].set(n_specset["species"]);


        BP_PLUGIN_INFO("Specset Info for "
                       << mesh_specset_name
                       << " : " << m_specset_info[mesh_specset_name].to_yaml())

        avtSpeciesMetaData *smd = new avtSpeciesMetaData(
            mesh_specset_name, // The name of the species
            mesh_topo_name,    // The name of the mesh the species is defined on.
            mesh_matset_name,  // The name of the material the species is defined on.
            nmat,              // The number of materials in the matset.
            num_species,       // The number of species for each material.
            species_names);    // The name of each species for each material.
        smd->validVariable = true;
        md->Add(smd);
    }
}


// ****************************************************************************
// helper method used to add expression meta data for a blueprint mesh.
//
// Mark C. Miller, Wed May  6 12:26:33 PDT 2020
// ****************************************************************************
static void
AddBlueprintExpressionMetadata(avtDatabaseMetaData *md, string const &mesh_name,
    const Node &n_mesh_info)
{
    if (!n_mesh_info.has_child("expressions"))
        return;

    BP_PLUGIN_INFO("adding expressions for " <<  mesh_name);

    NodeConstIterator exprs_itr = n_mesh_info["expressions"].children();

    while (exprs_itr.has_next())
    {
        const Node &n_expr = exprs_itr.next();

        if (n_expr.has_child("consumer") &&
            !StringHelpers::CaseInsensitiveEqual(n_expr["consumer"].as_string().c_str(), "visit"))
            continue;

        string expname = exprs_itr.name();
        string exp_topo_name = n_expr["topology"].as_string();
        string exp_mesh_name = mesh_name + "_" + exp_topo_name;
        string expname_wmesh = exp_mesh_name + "/" + expname;
        int ncomps = n_expr["number_of_components"].to_int();

        // Lookup what we already know about this mesh from database metadata
        const avtMeshMetaData *mmd = md->GetMesh(mesh_name + "_" + exp_topo_name);
        int ndims = mmd ? mmd->spatialDimension : 0;

        Expression expr;
        expr.SetName(expname_wmesh);
        if (ncomps == 1)
            expr.SetType(Expression::ScalarMeshVar);
        else if (ndims == 2 && ncomps == 2)
            expr.SetType(Expression::VectorMeshVar);
        else if (ndims == 2 && ncomps == 3)
            expr.SetType(Expression::SymmetricTensorMeshVar);
        else if (ndims == 2 && ncomps == 4)
            expr.SetType(Expression::TensorMeshVar);
        else if (ndims == 3 && ncomps == 3)
            expr.SetType(Expression::VectorMeshVar);
        else if (ndims == 3 && ncomps == 6)
            expr.SetType(Expression::SymmetricTensorMeshVar);
        else if (ndims == 3 && ncomps == 9)
            expr.SetType(Expression::TensorMeshVar);
        else
            expr.SetType(Expression::ArrayMeshVar);

        // Find all variable references in defn and replace with
        // parent mesh name. 

        std::string defn = n_expr["definition"].as_string();
        size_t ltidx=0;
        while ((ltidx = defn.find('<', ltidx)) != std::string::npos)
        {
            size_t gtidx = defn.find('>', ltidx+1);
            std::string vname(defn, ltidx+1, gtidx-ltidx-1);
            size_t slash = vname.find('/', 0);

            if (slash)
            {
                std::string b4slash(vname, 0, slash);

                // Skip this variable if it already has a meshname
                if (md->GetMesh(b4slash))
                {
                    ltidx += vname.length() + 1;
                    continue;
                }
            }

            defn.erase(ltidx+1, gtidx-ltidx-1);
            defn.insert(ltidx+1, exp_mesh_name + "/" + vname);
            ltidx += exp_mesh_name.length() + 1 + vname.length() + 1;
        }

        expr.SetDefinition(defn);
        md->AddExpression(&expr);
    }
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadRootFile
//
//  Purpose: Read contents of the root file
//
//  Programmer: cyrush
//  Creation:   Fri Dec  8 14:55:23 PST 2017
//
//  Modifications:
//    Cyrus Harrison, Fri Aug 24 14:01:50 PDT 2018
//    Add extra check for valid HDF5 file and allow plugin to be used if
//    any valid mesh index is found.
//
//    Cyrus Harrison, Mon Mar  9 15:15:55 PDT 2020
//    Refactor how index info is read to reduce I/O.
//
//    Cyrus Harrison, Tue Dec 13 12:17:17 PST 2022
//    Change index check to look for blueprint_index w/o fully reading it.
//
// ****************************************************************************
void
avtBlueprintFileFormat::ReadRootFile()
{
        //
        // Read root file using conduit::relay
        //
        // Note:
        // we only want to read the metadata portions
        // since folks can pack meshes into the root file as well :-)
        //

        string root_fname = GetFilename();

        BP_PLUGIN_INFO("Opening root file " << root_fname);

        int error = 0;

        // first figure out the protocol of the root file
        // this may be different than the protocol of the actual data files

        // assume hdf5, but check json file

        std::string root_protocol = "hdf5";
        std::string error_msg = "";

// only check on proc-0
#ifdef PARALLEL
        if (PAR_Rank() == 0)
#endif
        {
            conduit::relay::io::identify_file_type(root_fname,root_protocol);

            if(root_protocol.find("hdf5") != std::string::npos)
            {
               if(conduit::relay::io::is_hdf5_file(root_fname))
               {
                   // fast fail check for if this is a valid blueprint root file
                   try
                   {
                        relay::io::IOHandle root_hnd;
                        Node open_opts;
                        open_opts["mode"] = "r";
                        root_hnd.open(root_fname, root_protocol, open_opts);
                        if(  !root_hnd.is_open() || !root_hnd.has_path("blueprint_index"))
                        {
                            error_msg = "Root file: "
                                        + root_fname 
                                        + " does not contain `blueprint_index`.";
                            error = 1;
                        }
                        root_hnd.close();
                    }
                    catch(conduit::Error &e)
                    {
                        error_msg = e.message();
                        error = 1;
                    }
                }
                else
                {
                    error_msg = root_fname + " is not a valid HDF5 file.\n" +
                              " Cannot open with 'hdf5' protocol.";
                    error = 1;
                }
            }
        }

        if(root_protocol == "unknown")
        {
            error_msg = root_fname + " could not determine root file type (hdf5, json, or yaml )\n";
            error = 1;
        }

// check for error reading root file
#ifdef PARALLEL
        Node n_in, n_out;
        n_in.set(error);
        conduit::relay::mpi::sum_all_reduce(n_in,
                                            n_out,
                                            VISIT_MPI_COMM);

        error = n_out.to_int();
#endif
       if(error != 0)
       {
           BP_PLUGIN_EXCEPTION1(InvalidFilesException,
                                "Error reading root file: '" << root_fname<<"'. "
                                <<error_msg);
       }

#ifdef PARALLEL
        // only read on rank 0 and broadcast to everyone else
        if (PAR_Rank() == 0)
        {
            ReadRootIndexItems(root_fname,root_protocol,m_root_node);
        }

        conduit::relay::mpi::broadcast_using_schema(m_root_node,
                                                    0,
                                                    VISIT_MPI_COMM);
#else
        ReadRootIndexItems(root_fname,root_protocol,m_root_node);
#endif

        if(!m_root_node.has_child("blueprint_index"))
        {
            BP_PLUGIN_EXCEPTION1(InvalidFilesException,
                                 "Root file missing 'blueprint_index'");
        }

        NodeConstIterator itr = m_root_node["blueprint_index"].children();
        Node n_verify_info;

        bool any_index_ok = false;

        while(itr.has_next())
        {
            const Node &curr = itr.next();
            std::string mesh_name = itr.name();
            if( blueprint::mesh::index::verify(curr,
                                               n_verify_info[mesh_name]))
            {
                any_index_ok = true;

                debug5 << "Success: Mesh Blueprint Index verify for "
                      << mesh_name << std::endl;
            }
        }

        if(!any_index_ok)
        {
            BP_PLUGIN_EXCEPTION1(InvalidFilesException,
                                 "Failed to find a valid Mesh Blueprint Index\n"
                                 << n_verify_info.to_yaml());
        }
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::ReadRootIndexItems
//
//  Purpose: Read only the metadata portions of the root file that VisIt needs.
//
//  Programmer: cyrush
//  Creation:   Mon Mar  9 12:24:16 PDT 2020
//
//  Modifications:
//   Cyrus Harrison, Fri Jul 16 15:39:57 PDT 2021
//   Bugfix: When using relay::io::IOHandle, always open in read only mode
//
//   Cyrus Harrison, Tue Dec 13 12:04:01 PST 2022
//   Remove special case for HDF5 and use relay::io::IOHandle for all cases
//
//   Cyrus Harrison, Thu Nov  2 11:19:32 PDT 2023
//   Added call to AugmentBlueprintIndex to identify meshes with HO
//   volume fractions and auto a matset for them.
//
// ****************************************************************************
void
avtBlueprintFileFormat::ReadRootIndexItems(const std::string &root_fname,
                                           const std::string &root_protocol,
                                           conduit::Node &root_info)
{
    // list of names we want to read from the root file
    conduit::Node index_names;
    index_names.append() = "blueprint_index";
    index_names.append() = "file_pattern";
    index_names.append() = "tree_pattern";
    index_names.append() = "number_of_trees";
    index_names.append() = "number_of_files";
    index_names.append() = "protocol";

    // We don't want to return everything b/c the index meta data is
    // broadcasted to all ranks and is also printed in debug 5 logs
    // so we still filter what is pulled out here
    relay::io::IOHandle root_hnd;
    Node open_opts;
    open_opts["mode"] = "r";
    root_hnd.open(root_fname, root_protocol, open_opts);

    // loop over all names and copy them to the output node
    NodeConstIterator itr = index_names.children();
    while(itr.has_next())
    {
        std::string curr_idx_name = itr.next().as_string();
        if(root_hnd.has_path(curr_idx_name))
        {
            root_hnd.read(curr_idx_name,
                          root_info[curr_idx_name]);
        }
    }

    root_hnd.close();

    AugmentBlueprintIndex(root_info["blueprint_index"]);
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::AugmentBlueprintIndex
//
//  Purpose: Identify special case for high order materials and create a matset
//
//  Programmer: cyrush
//  Creation:   Wed Nov  1 15:10:15 PDT 2023
//
//
// ****************************************************************************
void
avtBlueprintFileFormat::AugmentBlueprintIndex(conduit::Node &blueprint_index)
{
    // loop over all meshes and see if we have high order volume fractions
    // that aren't published as a `matset`
    NodeIterator mesh_itr = blueprint_index.children();
    while(mesh_itr.has_next())
    {
        Node &mesh = mesh_itr.next();
        std::string mesh_name = mesh.name();
        index_t num_volfracs = 0;
        if(!mesh.has_path("matsets/matset") && mesh.has_child("fields"))
        {
            // walk index and look for the field volume_fraction_001
            if(mesh.has_path("fields/volume_fraction_001"))
            {
                BP_PLUGIN_INFO("Found HO material volume fractions for mesh: " << mesh_name);
                // we have materials !!!!
                // count to find number of vol fracs
                NodeConstIterator fields_itr = mesh["fields"].children();
                while(fields_itr.has_next())
                {
                    const Node &field = fields_itr.next();
                    if(field.name().find("volume_fraction_") !=std::string::npos)
                    {
                        num_volfracs++;
                    }
                }
            }

            if(num_volfracs > 0)
            {
                BP_PLUGIN_INFO("Matset for mesh " << mesh_name  <<
                                "has " << num_volfracs << " material regions.");
                // create matset
                Node &mset = mesh["matsets/matset"];
                mset["topology"] = mesh["fields/volume_fraction_001/topology"];
                // fake it
                mset["path"] = mesh["fields/volume_fraction_001/path"];
                Node &mats_list = mset["materials"];
                for(index_t idx=1;idx<num_volfracs+1;idx++)
                {
                    mats_list[conduit_fmt::format("{:03d}",idx)] = idx;
                }

                BP_PLUGIN_INFO("Auto generated matset for mesh: " << mesh_name <<
                               std::endl << mesh["matsets/matset"].to_yaml());
            }
        }
    }
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
//
//  Modifications:
//    Cyrus Harrison, Mon Mar  9 15:45:17 PDT 2020
//    Use explicit map from registered mesh name to bp mesh and topo names.
//
//    Cyrus Harrison, Tue Dec 13 12:17:17 PST 2022
//    Refactor to pass mesh name to tree cache for part map style index
//    support.
//
//    Eric Brugger, Fri Feb 24 14:57:15 PST 2023
//    Added a stop timer to a catch block to avoid the case of not calling
//    stop timer when an exception occurs.
// 
//    Justin Privitera, Wed Mar 22 16:09:52 PDT 2023
//    Bookkeeping for 1D curves.
// 
//    Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//    Fixed warnings.
// 
//    Justin Privitera, Fri Sep 27 11:51:59 PDT 2024
//    Added support for species.
//
// ****************************************************************************

void
avtBlueprintFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *)
{
    BP_PLUGIN_INFO("Begin avtBlueprintFileFormat::PopulateDatabaseMetaData");

    int t_pop_md = visitTimer->StartTimer();

    // clear any mfem mesh mappings
    m_mfem_mesh_map.clear();
    // clear any 1D curve names
    m_curve_names.clear();
    // clear full mesh to bp mesh and topo name map
    m_mesh_and_topo_info.reset();
    // clear full matset to info map
    m_matset_info.reset();
    // clear full specset to info map
    m_specset_info.reset();

    try
    {
        ReadRootFile();

        BP_PLUGIN_INFO("Root file contents" << endl << m_root_node.to_yaml());

        m_protocol = "hdf5";

        if(m_root_node.has_child("protocol"))
        {
            m_protocol = m_root_node["protocol/name"].as_string();
        }

        if (DebugStream::Level5())
        {
            m_root_node.to_yaml_stream(DebugStream::Stream5());
        }

        BP_PLUGIN_INFO("Using protocol: " << m_protocol);

        m_tree_cache->SetProtocol(m_protocol);

        string root_fname = GetFilename();
        string root_dir = FileFunctions::Dirname(root_fname);

        m_tree_cache->SetRootDir(root_dir);

        // simple / legacy index case
        if(m_root_node.has_child("file_pattern"))
        {

            string file_pattern = m_root_node["file_pattern"].as_string();

            // if file_pattern isn't an abs path, it needs to be relative to the
            // the root file
            if(file_pattern[0] !=  VISIT_SLASH_STRING[0])
            {
                file_pattern  = root_dir + string(VISIT_SLASH_STRING) + file_pattern;
            }

            m_tree_cache->SetFilePattern(file_pattern);
            m_tree_cache->SetTreePattern(m_root_node["tree_pattern"].as_string());
            m_tree_cache->SetNumberOfFiles(m_root_node["number_of_files"].to_int());
            m_tree_cache->SetNumberOfTrees(m_root_node["number_of_trees"].to_int());
        }

        NodeConstIterator itr = m_root_node["blueprint_index"].children();

        while (itr.has_next())
        {
            const Node &n = itr.next();
            // check for per mesh part map case
            if(n.has_path("state/partition_pattern"))
            {
                BP_PLUGIN_INFO("Adding Blueprint Index Partition Map for "
                               << itr.name());

                if(n.has_path("state/partition_map"))
                {
                    m_tree_cache->AddMeshParitionMap(itr.name(),
                                                     n["state/partition_pattern"].as_string(),
                                                     n["state/partition_map"]);
                }
                else // special case (single file)
                {
                    m_tree_cache->AddMeshParitionMap(itr.name(),
                                                     n["state/partition_pattern"].as_string());
                }
            }

            AddBlueprintMeshAndFieldMetadata(metadata, itr.name(), n);
            AddBlueprintMaterialsMetadata(metadata, itr.name(), n);
            AddBlueprintSpeciesMetadata(metadata, itr.name(), n);
        }

        // Process all expressions *after* all fields. This
        // helps filter expression content relative to all known meshes.
        itr = m_root_node["blueprint_index"].children();

        while (itr.has_next())
        {
            const Node &n = itr.next();
            AddBlueprintExpressionMetadata(metadata, itr.name(), n);
        }
    }
    catch(conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "Populate Database MetaData: " << endl
                    << e.message();
        visitTimer->StopTimer(t_pop_md,"PopulateDatabaseMetaData");
        EXCEPTION1(InvalidFilesException, err_oss.str());
    }


    visitTimer->StopTimer(t_pop_md,"PopulateDatabaseMetaData");
    BP_PLUGIN_INFO("End avtBlueprintFileFormat::PopulateDatabaseMetaData");
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetCycle
//
//  Purpose:
//      Returns if we have the current cycle.
//
//  Programmer: Cyrus Harrison
//  Creation:   Thu Aug 11 10:15:05 PDT 2016
// 
//  Modifications:
//     Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//     Encased in try-catch block.
//
// ****************************************************************************
int
avtBlueprintFileFormat::GetCycle()
{
    try
    {
        // VisIt doesn't support diff times / cycles for meshes in STMD
        // we loop over all meshes and return the first valid cycle
        NodeConstIterator itr = m_root_node["blueprint_index"].children();

        while(itr.has_next())
        {
            const Node &mesh = itr.next();
            if(mesh.has_path("state/cycle"))
            {
                return mesh["state/cycle"].to_int();
            }
        }

        return avtFileFormat::INVALID_CYCLE;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetCycle: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetTime
//
//  Purpose:
//      Returns if we have the current cycle.
//
//  Programmer: Cyrus Harrison
//  Creation:   Thu Aug 11 10:15:05 PDT 2016
// 
//  Modifications:
//     Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//     Encased in try-catch block.
//
// ****************************************************************************
double
avtBlueprintFileFormat::GetTime()
{
    try
    {
        // VisIt doesn't support diff times / cycles for meshes in STMD
        // we loop over all meshes and return the first valid time

        NodeConstIterator itr = m_root_node["blueprint_index"].children();

        while(itr.has_next())
        {
            const Node &mesh = itr.next();
            if(mesh.has_path("state/time"))
            {
                return mesh["state/time"].to_double();
            }
        }

        return avtFileFormat::INVALID_TIME;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetTime: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
//
//  Modifications:
//    Cyrus Harrison, Tue Mar 10 13:15:26 PDT 2020
//    Support empty mesh case.
//
//    Justin Privitera, Wed Mar 23 12:24:49 PDT 2022
//    Passes domain to MeshToVTK.
// 
//    Justin Privitera, Mon Apr 11 18:20:19 PDT 2022
//    `m_new_refine` is passed to `RefineMeshToVTK` so it can choose the 
//    appropriate MFEM LOR scheme.
// 
//    Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//    Encased in try-catch block.
//
//    Justin Privitera, Wed Mar 22 16:09:52 PDT 2023
//    Added support for 1D curves.
//
//    Brad Whitlock, Mon May 22 17:29:05 PDT 2023
//    I added some code to clear the mesh's material from the cache if it
//    has an obsolete LOD.
// 
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
//
//    Justin Privitera, Fri May  3 09:55:25 PDT 2024
//    Change how we fetch ndims so that we do not add "values" leaf to 
//    uniform coordsets.
// 
// ****************************************************************************

vtkDataSet *
avtBlueprintFileFormat::GetMesh(int domain, const char *abs_meshname)
{
    try
    {
        BP_PLUGIN_INFO("GetMesh: "
                       << abs_meshname << " [domain " << domain << "]");

        // read mesh data into conduit tree
        Node data;
        string abs_meshname_str(abs_meshname);
        string abs_1d_curve_field_name{""};

        bool curve_1d = m_curve_names.find(abs_meshname_str) != m_curve_names.end();
        if (curve_1d)
        {
            // This means that abs_meshname is of the form:
            // "mesh_name/field_name", so we need to separate it
            abs_meshname_str = abs_meshname_str.substr(0, abs_meshname_str.find('/'));
            abs_1d_curve_field_name = FileFunctions::Basename(abs_meshname);
        }

        // reads a single mesh into a blueprint conforming output
        ReadBlueprintMesh(domain, abs_meshname_str, data);
        // check for empty mesh case
        if(data.dtype().is_empty())
        {
            // support empty mesh case by returning NULL
            return NULL;
        }

        Node verify_info;
        if(!blueprint::mesh::verify(data, verify_info))
        {
            BP_PLUGIN_INFO("blueprint::mesh::verify failed for mesh "
                           << abs_meshname << " [domain " << domain << "]" << endl
                           << "Verify Info " << endl
                           << verify_info.to_yaml() << endl
                           << "Data Schema " << endl
                           << data.schema().to_yaml());

            BP_PLUGIN_INFO("warning: "
                           "avtBlueprintFileFormat::GetMesh returning NULL "
                           << abs_meshname
                           << " [domain " << domain << "]"
                           << " will be missing" << endl);
            // TODO: Should we throw an error instead of blanking the domain?
            return NULL;
        }

        BP_PLUGIN_INFO(data.schema().to_yaml());

        // prepare result vtk dataset
        vtkDataSet *res = NULL;


        string mesh_name;
        string topo_name;

        FetchMeshAndTopoNames(abs_meshname_str,
                              mesh_name,
                              topo_name);

        BP_PLUGIN_INFO("mesh name: " << mesh_name);
        BP_PLUGIN_INFO("topo name: " << topo_name);

        const conduit::Node &n_coords = data["coordsets"][0];
        int ndims;
        if (n_coords.has_child("values"))
        {
            ndims = static_cast<int>(n_coords["values"].number_of_children());
        }
        else
        {
            // if you don't have values, you are a uniform coordset so you have dims
            ndims = static_cast<int>(n_coords["dims"].number_of_children());
        }

        // check for the mfem case
        if( m_mfem_mesh_map[topo_name] )
        {
            BP_PLUGIN_INFO("mesh  " << topo_name << " is a mfem mesh");

            // If there are any materials associated with this mesh in the cache
            // that do not match the current refinement level, purge them so the
            // database will re-read them at the selected level of detail.
            auto it = m_mfem_material_map.find(topo_name);
            if(it != m_mfem_material_map.end())
            {
                if(it->second.second != (m_selected_lod + 1))
                {
                    BP_PLUGIN_INFO("Material " << it->second.first << " has LOD "
                        << it->second.second << ". It needs to be purged.");
                    cache->ClearVariablesWithString(it->second.first);
                    m_mfem_material_map.erase(it);
                }
            }

            // use mfem to refine and create a vtk dataset
            mfem::Mesh *mesh = avtConduitBlueprintDataAdaptor::BlueprintToMFEM::MeshToMFEM(data);
            res = avtMFEMDataAdaptor::RefineMeshToVTK(mesh,
                                                      domain,
                                                      m_selected_lod+1,
                                                      m_new_refine);

            // cleanup the mfem mesh
            delete mesh;
        }
        else if (ndims == 1)
        {
            BP_PLUGIN_INFO("mesh  " << topo_name << " is a 1D curve");

            Node field_1d;
            ReadBlueprintField(domain, abs_meshname, field_1d);
            res = avtConduitBlueprintDataAdaptor::BlueprintToVTK::Curve1DToVTK(n_coords, field_1d);
        }
        else
        {
            BP_PLUGIN_INFO("mesh  " << topo_name << " is a standard mesh");
            // construct a vtk dataset directly from blueprint data
            // in a conduit tree
            res = avtConduitBlueprintDataAdaptor::BlueprintToVTK::MeshToVTK(domain, data);
        }

        BP_PLUGIN_INFO("avtBlueprintFileFormat::GetMesh Done");
        return res;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetMesh: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
//
//  Modifications:
//    Cyrus Harrison, Wed Mar 11 10:42:22 PDT 2020
//    Allow empty domains.
// 
//    Justin Privitera, Wed Mar 23 12:21:50 PDT 2022
//    Now handles polyhedral/polygonal meshes and converts them
//    to tetrahedral/triangular meshes using conduit.
//
//    Cyrus Harrison, Mon Mar 28 12:14:20 PDT 2022
//    Use conduit version check for polytopal support.
//
//    Justin Privitera, Mon May 23 20:19:52 PDT 2022
//    `m_new_refine` is passed to `RefineGridFunctionToVTK` so it can choose
//    the appropriate MFEM LOR scheme.
// 
//    Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//    Encased in try-catch block.
// 
//    Justin Privitera, Tue Aug 23 14:40:24 PDT 2022
//    Removed `CONDUIT_HAVE_PARTITION_FLATTEN` check.
// 
//    Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//    Fixed warnings.
//
//    Brad Whitlock, Tue Dec 19 14:45:13 PST 2023
//    Add MFEM refinement for certain fields that are low-order but live on a
//    mesh that has been refined. This prevents the field from having too few
//    values for the refined mesh.
// 
//    Justin Privitera, Sat Jun 29 14:22:21 PDT 2024
//    Handle the polytopal mixed topology case.
//
// ****************************************************************************

vtkDataArray *
avtBlueprintFileFormat::GetVar(int domain, const char *abs_varname)
{
    try
    {
        BP_PLUGIN_INFO("GetVar: " << abs_varname << " [domain " << domain << "]");

        vtkDataArray *res = NULL;

        string abs_varname_str(abs_varname);

        // check for special var "element_coloring"
        if(abs_varname_str.find("element_coloring") != std::string::npos)
        {
            // TODO: we currently have to replace colons, brackets, etc
            // to fetch from the mesh metadata, is there a standard helper
            // util in VisIt that can take care of this for us?
            string abs_meshname = metadata->MeshForVar(sanitize_var_name(abs_varname_str));

            // element coloring is generated from the mesh
            Node n_mesh;
            // read the mesh data
            ReadBlueprintMesh(domain, abs_meshname, n_mesh);
            // check for empty mesh case
            if(n_mesh.dtype().is_empty())
            {
                // support empty mesh case by returning NULL
                return NULL;
            }

            Node verify_info;
            if(!blueprint::mesh::verify(n_mesh,verify_info))
            {
                BP_PLUGIN_INFO("blueprint::mesh::verify failed for field "
                               << abs_meshname << " [domain " << domain << "]" << endl
                               << "Verify Info " << endl
                               << verify_info.to_yaml() << endl
                               << "Data Schema " << endl
                               << n_mesh.schema().to_yaml());
                return NULL;
            }

            // create an mfem mesh
            mfem::Mesh *mesh = avtConduitBlueprintDataAdaptor::BlueprintToMFEM::MeshToMFEM(n_mesh);
            // refine the coloring to a vtk data array
            res = avtMFEMDataAdaptor::RefineElementColoringToVTK(mesh,
                                                                 domain,
                                                                 m_selected_lod+1);
            // clean up the mfem mesh
            delete mesh;

            // return the coloring result
            return res;
        }

        // else, normal field case

        Node n_field;
        Node intermediate_data;
        Node *field_ptr = &n_field;

        ReadBlueprintField(domain,abs_varname_str,*field_ptr);

        // check for empty field case
        if(field_ptr->dtype().is_empty())
        {
            // support empty field case by returning NULL
            return NULL;
        }

        Node field_verify_info;
        if(!blueprint::mesh::field::verify(*field_ptr,field_verify_info))
        {
            BP_PLUGIN_INFO("blueprint::mesh::field::verify failed for field "
                           << abs_varname_str << " [domain " << domain << "]" << endl
                           << "Verify Info " << endl
                           << field_verify_info.to_yaml() << endl
                           << "Data Schema " << endl
                           << field_ptr->schema().to_yaml());
            return NULL;
        }

        // if we have an association, this is a standard field
        if(field_ptr->has_child("association"))
        {
            string abs_meshname = metadata->MeshForVar(sanitize_var_name(abs_varname_str));
            // element coloring is generated from the mesh
            Node n_mesh;
            // read the mesh data
            ReadBlueprintMesh(domain, abs_meshname, n_mesh);
            // check for empty mesh case
            if(n_mesh.dtype().is_empty())
            {
                // support empty mesh case by returning NULL
                return NULL;
            }

            Node verify_info;
            if(!blueprint::mesh::verify(n_mesh,verify_info))
            {
                BP_PLUGIN_INFO("blueprint::mesh::verify failed for field "
                               << abs_meshname << " [domain " << domain << "]" << endl
                               << "Verify Info " << endl
                               << verify_info.to_yaml() << endl
                               << "Data Schema " << endl
                               << n_mesh.schema().to_yaml());
                return NULL;
            }

            string mesh_name, topo_name;

            FetchMeshAndTopoNames(abs_meshname,
                                  mesh_name,
                                  topo_name);

            BP_PLUGIN_INFO("mesh name: " << mesh_name);
            BP_PLUGIN_INFO("topo name: " << topo_name);

            Node &n_topo = n_mesh["topologies/" + topo_name];
            Node *topo_ptr = &n_topo;
            const Node &n_coords = n_mesh["coordsets/" + n_topo["coordset"].as_string()];
            int tdims = -1;
            if (n_topo.has_path("elements/shape"))
            {
                // Get the shape's topological dimension.
                conduit::blueprint::mesh::utils::ShapeType shape(n_topo);
                tdims = shape.dim;

                if (shape.is_polygonal() || shape.is_polyhedral())
                {
                    string field_name, coords_name;
                    // get name of coordset from topology
                    coords_name = n_topo["coordset"].as_string();

                    Node s2dmap, d2smap, options;
                    Node &side_coords = intermediate_data["side_mesh"]["coordsets/" + coords_name];
                    Node &side_topo = intermediate_data["side_mesh"]["topologies/" + topo_name];
                    Node &side_fields = intermediate_data["side_mesh"]["fields"];

                    field_name = FileFunctions::Basename(abs_varname_str);
                    n_mesh["fields/" + field_name].set_external(*field_ptr);

                    blueprint::mesh::topology::unstructured::generate_sides(
                      n_mesh["topologies/" + topo_name],
                      side_topo,
                      side_coords,
                      side_fields,
                      s2dmap,
                      d2smap,
                      options);

                    field_ptr = side_fields.fetch_ptr(field_name);
                    topo_ptr = intermediate_data.fetch_ptr("side_mesh/topologies/" + topo_name);
                }
                // polytopal mixed case
                else if (n_topo["elements/shape"].as_string() == "mixed" &&
                         (n_topo.has_child("subelements") || 
                          n_topo.has_path("elements/shape_map/polygonal")))
                {
                    const string field_name = FileFunctions::Basename(abs_varname_str);
                    // either we are making a polyhedral mesh or a polygonal mesh
                    const bool mesh_is_polyhedral = n_topo.has_child("subelements");

                    if (mesh_is_polyhedral && n_topo.has_path("elements/shape_map/polygonal"))
                    {
                        BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                             "The mixed polygonal and polyhedral mesh "
                                             "case is currently unsupported.");
                    }

                    // 
                    // step 1: threshold out the polytopal elements, placing 
                    // them in their own topology
                    // 
                    const std::string coordset_name = n_topo["coordset"].as_string();
                    Node &polytopal_mesh = intermediate_data["mixed_transformation/polytopal_mesh"];
                    
                    const int num_polytopal_elements = 
                        avtConduitBlueprintDataAdaptor::BlueprintToVTK::CreatePolytopalMeshFromMixedMesh(
                            n_coords, 
                            n_topo,
                            polytopal_mesh);

                    const Node &polytopal_topo = polytopal_mesh["topologies"][n_topo.name()];

                    int_accessor n_shapes = n_topo["elements"]["shapes"].value();

                    // we must create the field on the polytopal mesh ourselves
                    Node &polytopal_field = polytopal_mesh["fields"][field_name];
                    const bool elem_assoc = field_ptr->fetch("association").as_string() == "element";
                    auto field_entry_itr = field_ptr->children();
                    while (field_entry_itr.has_next())
                    {
                        const Node &field_entry = field_entry_itr.next();
                        const std::string entry_name = field_entry_itr.name();

                        // modify field values if we are element associated
                        // grab only the ones we need
                        if (entry_name == "values" && elem_assoc)
                        {
                            double_accessor n_field_values = n_field["values"].value();

                            polytopal_field["values"].set(DataType::float64(num_polytopal_elements));
                            double_array new_field_values = polytopal_field["values"].value();
                            int new_field_index = 0;

                            if (mesh_is_polyhedral)
                            {
                                for (int zoneid = 0; zoneid < n_shapes.dtype().number_of_elements(); zoneid ++)
                                {
                                    if (n_shapes[zoneid] == VTK_POLYHEDRON)
                                    {
                                        new_field_values[new_field_index] = n_field_values[zoneid];
                                        new_field_index ++;
                                    }
                                }
                            }
                            else
                            {
                                for (int zoneid = 0; zoneid < n_shapes.dtype().number_of_elements(); zoneid ++)
                                {
                                    if (n_shapes[zoneid] == VTK_POLYGON)
                                    {
                                        new_field_values[new_field_index] = n_field_values[zoneid];
                                        new_field_index ++;
                                    }
                                }
                            }
                        }
                        else
                        {
                            polytopal_field[entry_name].set_external(field_entry);
                        }
                    }

                    // 
                    // step 2: run the polytopal mesh through generate_sides
                    // 
                    Node &side_mesh = intermediate_data["mixed_transformation/side_mesh"];
                    Node &side_topo = side_mesh[n_topo.name()];
                    Node &side_coords = side_mesh[coordset_name];
                    Node &side_fields = side_mesh["fields"];
                    Node &s2dmap = intermediate_data["mixed_transformation/maps/s2dmap"];
                    Node &d2smap = intermediate_data["mixed_transformation/maps/d2smap"];
                    blueprint::mesh::topology::unstructured::generate_sides(
                        polytopal_topo,
                        side_topo,
                        side_coords,
                        side_fields,
                        s2dmap,
                        d2smap);

                    // 
                    // step 3: stitch the topology back together to create a 
                    // brand new mixed topology
                    // 

                    Node &new_mixed_topo = intermediate_data["mixed_transformation/new_mixed_topo"];

                    if (mesh_is_polyhedral)
                    {
                        if (side_topo["elements/shape"].as_string() != "tet")
                        {
                            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                                 "Generated elements for mixed polyhedral "
                                                 "topology must be tetrahedrons.");
                        }
                    }
                    else
                    {
                        if (side_topo["elements/shape"].as_string() != "tri")
                        {
                            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                                 "Generated elements for mixed polygonal "
                                                 "topology must be triangles.");
                        }
                    }

                    avtConduitBlueprintDataAdaptor::BlueprintToVTK::CreateMixedMeshFromSideAndMixedMeshes(
                        n_topo,
                        side_topo,
                        new_mixed_topo);

                    // 
                    // step 4: stitch the field back together
                    // 

                    Node &new_field = intermediate_data["mixed_transformation/fields/" + field_name];
                    field_entry_itr = side_fields[field_name].children();
                    while (field_entry_itr.has_next())
                    {
                        const Node &field_entry = field_entry_itr.next();
                        const std::string entry_name = field_entry_itr.name();

                        // modify field values if we are element associated
                        // grab only the ones we need
                        if (entry_name == "values" && elem_assoc)
                        {
                            double_accessor n_field_values = n_field["values"].value();
                            double_accessor side_field_values = side_fields[field_name]["values"].value();

                            // the number of elements is the number of non-polytopal elements plus the number 
                            // of new triangular/tetrahedral elements.
                            const int final_num_elems = n_shapes.number_of_elements() - num_polytopal_elements
                                                      + side_field_values.number_of_elements();

                            new_field["values"].set(DataType::float64(final_num_elems));
                            double_array new_field_values = new_field["values"].value();
                            int new_field_index = 0;

                            // grab original field values
                            if (mesh_is_polyhedral)
                            {
                                for (int zoneid = 0; zoneid < n_shapes.dtype().number_of_elements(); zoneid ++)
                                {
                                    // grab all the original shapes
                                    if (n_shapes[zoneid] != VTK_POLYHEDRON)
                                    {
                                        new_field_values[new_field_index] = n_field_values[zoneid];
                                        new_field_index ++;
                                    }
                                }
                            }
                            else
                            {
                                for (int zoneid = 0; zoneid < n_shapes.dtype().number_of_elements(); zoneid ++)
                                {
                                    // grab all the original shapes
                                    if (n_shapes[zoneid] != VTK_POLYGON)
                                    {
                                        new_field_values[new_field_index] = n_field_values[zoneid];
                                        new_field_index ++;
                                    }
                                }
                            }
                            
                            // get new field values
                            for (int zoneid = 0; zoneid < side_field_values.dtype().number_of_elements(); zoneid ++)
                            {
                                // grab all the new shapes
                                new_field_values[new_field_index] = side_field_values[zoneid];
                                new_field_index ++;
                            }
                        }
                        else
                        {
                            new_field[entry_name].set_external(field_entry);
                        }
                    }

                    // 
                    // step 5: update the reference to the field and topology to
                    // point to the new ones we created
                    // 
                    field_ptr = intermediate_data.fetch_ptr("mixed_transformation/fields/" + field_name);
                    topo_ptr = intermediate_data.fetch_ptr("mixed_transformation/new_mixed_topo");
                }
            }

            // Check whether the data need to be refined anyway.
            const bool meshIsHO = n_mesh.has_path("topologies/" + topo_name + "/grid_function");
            if(meshIsHO && (tdims >= 2) && ((m_selected_lod+1) > 1))
            {
                // The field data are low-order but need to be refined anyway since
                // the mesh got refined.

                // Wrap as HO-compatible field.
                conduit::Node fieldHO;
                fieldHO["topology"] = topo_name;
                if(field_ptr->fetch_existing("association").as_string() == "element")
                    fieldHO["basis"] = (tdims == 2) ? "L2_T2_2D_P0" : "L2_T2_3D_P0";
                else
                    fieldHO["basis"] = (tdims == 2) ? "H1_2D_P1" : "H1_3D_P1";
                fieldHO["values"].set_external(field_ptr->fetch_existing("values"));

                // create an mfem mesh
                mfem::Mesh *mesh = avtConduitBlueprintDataAdaptor::BlueprintToMFEM::MeshToMFEM(n_mesh);

                // create the grid fuction
                mfem::GridFunction *gf =  avtConduitBlueprintDataAdaptor::BlueprintToMFEM::FieldToMFEM(mesh,
                                                                                                       fieldHO);
                // refine the grid function into a vtk data array
                res = avtMFEMDataAdaptor::RefineGridFunctionToVTK(mesh,
                                                                  gf,
                                                                  m_selected_lod+1,
                                                                  m_new_refine);

                // cleanup mfem data
                delete gf;
                delete mesh;
            }
            else
            {
                // low-order case, use vtk
                res = avtConduitBlueprintDataAdaptor::BlueprintToVTK::FieldToVTK(*topo_ptr,
                                                                                 *field_ptr);
            }
        }
        // if we have a basis, this field is actually an mfem grid function
        else if(field_ptr->has_child("basis"))
        {
            // TODO: we currently have to replace colons, brackets, etc
            // to fetch from the mesh metadata, is there a standard helper
            // util in VisIt that can take care of this for us?
            string abs_meshname = metadata->MeshForVar(sanitize_var_name(abs_varname_str));

            // the grid function needs the mesh in order to refine

            // read the mesh data
            Node n_mesh;
            ReadBlueprintMesh(domain, abs_meshname, n_mesh);
            // check for empty mesh case
            if(n_mesh.dtype().is_empty())
            {
                // support empty mesh case by returning NULL
                return NULL;
            }

            Node verify_info;
            if(!blueprint::mesh::verify(n_mesh,verify_info))
            {
                BP_PLUGIN_INFO("blueprint::mesh::verify failed for field "
                               << abs_meshname << " [domain " << domain << "]" << endl
                               << "Verify Info " << endl
                               << verify_info.to_yaml() << endl
                               << "Data Schema " << endl
                               << n_mesh.schema().to_yaml());
                return NULL;
            }

            // create an mfem mesh
            mfem::Mesh *mesh = avtConduitBlueprintDataAdaptor::BlueprintToMFEM::MeshToMFEM(n_mesh);

            // create the grid fuction
            mfem::GridFunction *gf =  avtConduitBlueprintDataAdaptor::BlueprintToMFEM::FieldToMFEM(mesh,
                                                                                 *field_ptr);
            // refine the grid function into a vtk data array
            res = avtMFEMDataAdaptor::RefineGridFunctionToVTK(mesh,
                                                              gf,
                                                              m_selected_lod+1,
                                                              m_new_refine);

            // cleanup mfem data
            delete gf;
            delete mesh;
        }

        // check to see if we have matset_values multi-material data
        // if so, provide them as AUXILIARY_DATA_MIXED_VARIABLE
        if(field_ptr->has_child("matset"))
        {
            // get mesh and topo for var:
            //  - abs_varname_str
            
            // replace colons, etc
            abs_varname_str = sanitize_var_name(abs_varname_str);
            string abs_meshname = metadata->MeshForVar(abs_varname_str);

            BP_PLUGIN_INFO("field " << abs_varname << " is defined on mesh " << abs_meshname);

            string mesh_name;
            string topo_name;
            FetchMeshAndTopoNames(abs_meshname,
                                  mesh_name,
                                  topo_name);

            std::string matset_name = (*field_ptr)["matset"].as_string();
            std::string mat_name = mesh_name + "_" + topo_name + "_" + matset_name;

            BP_PLUGIN_INFO("mesh name: " << mesh_name);
            BP_PLUGIN_INFO("topo name: " << topo_name);
            BP_PLUGIN_INFO("materials name: " << mat_name);

            // get the matset
            // do the xform on the matset and the matset_values
            Node n_matset;
            ReadBlueprintMatset(domain,
                                mat_name,
                                n_matset);

            Node n_silo_matset;
            conduit::blueprint::mesh::field::to_silo(*field_ptr,
                                                     n_matset,
                                                     n_silo_matset);

            int mix_len  = static_cast<int>(n_silo_matset["field_mixvar_values"].dtype().number_of_elements());

            float *mixvals_ptr = NULL;
            if(n_silo_matset["field_mixvar_values"].dtype().is_float())
            {
                mixvals_ptr = n_silo_matset["field_mixvar_values"].as_float_ptr();
            }
            else
            {
                n_silo_matset["field_mixvar_values"].to_float_array(
                                            n_silo_matset["field_mixvar_values_float"]);
                mixvals_ptr = n_silo_matset["field_mixvar_values_float"].as_float_ptr();
            }

            avtMixedVariable *mvar = new avtMixedVariable(mixvals_ptr,
                                                          mix_len,
                                                          abs_varname_str);
            void_ref_ptr mvar_ref = void_ref_ptr(mvar, avtMixedVariable::Destruct);
            cache->CacheVoidRef(abs_varname_str.c_str(),
                                AUXILIARY_DATA_MIXED_VARIABLE,
                                timestep,
                                domain,
                                mvar_ref);
        }

        return res;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetVar: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }    
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: harrison37 -- generated by xml2avt
//  Creation:   Wed Jun 15 16:25:28 PST 2016
//
// ****************************************************************************

vtkDataArray *
avtBlueprintFileFormat::GetVectorVar(int domain, const char *varname)
{
    // vector vars can simply use the normal GetVar logic
    return GetVar(domain,varname);
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data from a Blueprint Database.
//
//  Arguments:
//      var        The variable of interest.
//      domain     The domain of interest.
//      type       The type of auxiliary data.
//      <unnamed>  The arguments for that -- not used
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: Cyrus Harrison
//  Creation:   December 8, 2020
//
//  Modifications:
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
// 
//     Justin Privitera, Fri Sep 27 11:51:59 PDT 2024
//     Added support for species.
//
// ****************************************************************************
void *
avtBlueprintFileFormat::GetAuxiliaryData(const char *var,
                                         int domain,
                                         const char *type,
                                         void * /* args (unused) */,
                                         DestructorFunction &df)
{
    void *rv = NULL;

    if (strcmp(type, AUXILIARY_DATA_MATERIAL) == 0)
    {
        rv = static_cast<void *>(GetMaterial(domain, var));
        df = avtMaterial::Destruct;
    }

    if (strcmp(type, AUXILIARY_DATA_SPECIES) == 0)
    {
        rv = static_cast<void *>(GetSpecies(domain, var));
        df = avtSpecies::Destruct;
    }

    return rv;
}


// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetMaterial
//
//  Purpose:
//      Gets the matset object from the Blueprint database and converts it 
//      into an avtMaterial.
//
//  Arguments:
//      domain     The domain of interest.
//      mat_name   The material of interest.
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: Cyrus Harrison
//  Creation:   December 8, 2020
//
//  Modifications:
//     Justin Privitera, Wed Aug 24 11:08:51 PDT 2022
//     Encased in try-catch block.
//
//     Justin Privitera, Tue Sep 19 11:36:45 PDT 2023
//     Get the material numbers and use a different avtMaterial constructor.
// 
//     Justin Privitera, Thu Oct 26 12:26:32 PDT 2023
//     Fixed warnings.
// 
//     Justin Privitera, Mon Feb  5 14:14:19 PST 2024
//     Removed unnecessary material numbers logic now that we have a new 
//     Conduit.
// 
//     Justin Privitera, Wed Feb 14 11:37:06 PST 2024
//     Present material ids alongside material names.
// 
//     Justin Privitera, Fri Mar 15 15:56:13 PDT 2024
//     Revert previous change.
//
// ****************************************************************************
avtMaterial *
avtBlueprintFileFormat::GetMaterial(int domain,
                                    const char *mat_name)
{
    try
    {
        BP_PLUGIN_INFO("avtBlueprintFileFormat::GetMaterial " 
                        << domain << " "
                        << mat_name);

        Node n_matset;
        ReadBlueprintMatset(domain,
                            mat_name,
                            n_matset);

        std::vector<std::string> matnames = n_matset["matnames"].child_names();
        // package up char ptrs
        std::vector<const char *> matnames_ptrs;
        for (const auto &matname : matnames)
            matnames_ptrs.push_back(matname.c_str());
        auto names = const_cast<char **>(matnames_ptrs.data());

        // use to_silo util to convert from bp to the mixslot rep
        // that silo and visit use

        Node n_silo_matset;
        conduit::blueprint::mesh::matset::to_silo(n_matset,
                                                  n_silo_matset);

        int nmats = static_cast<int>(matnames.size());
        int nzones = static_cast<int>(n_silo_matset["matlist"].dtype().number_of_elements());
        int *matlist  = NULL;
        int *mix_mat  = NULL;
        int *mix_next = NULL;

        // get the material numbers
        std::vector<int> matnos;
        auto matmap_itr = n_silo_matset["material_map"].children();
        while (matmap_itr.has_next())
        {
            const Node &n_mat = matmap_itr.next();
            matnos.push_back(n_mat.to_int());
        }

        // we need int ptrs for the avtMaterial object,
        // convert if needed

        Node n_tmp;
        if(!n_silo_matset["matlist"].dtype().is_int())
        {
            n_silo_matset["matlist"].to_int_array(n_tmp["matlist"]);
            n_silo_matset["mix_mat"].to_int_array(n_tmp["mix_mat"]);
            n_silo_matset["mix_next"].to_int_array(n_tmp["mix_next"]);
            matlist  = n_tmp["matlist"].as_int_ptr();
            mix_mat  = n_tmp["mix_mat"].as_int_ptr();
            mix_next = n_tmp["mix_next"].as_int_ptr();
        }
        else
        {
            matlist  = n_silo_matset["matlist"].as_int_ptr();
            mix_mat  = n_silo_matset["mix_mat"].as_int_ptr();
            mix_next = n_silo_matset["mix_next"].as_int_ptr();
        }

        int mix_len  = static_cast<int>(n_silo_matset["mix_mat"].dtype().number_of_elements());

        float *mix_vf = NULL;
        if(n_silo_matset["mix_vf"].dtype().is_float())
        {
            mix_vf = n_silo_matset["mix_vf"].as_float_ptr();
        }
        else
        {
            n_silo_matset["mix_vf"].to_float_array(n_silo_matset["mix_vf_float"]);
            mix_vf = n_silo_matset["mix_vf_float"].as_float_ptr();
        }

        const std::string domain_name = std::to_string(domain);

        avtMaterial *mat = new avtMaterial(nmats,               // The number of materials
                                           matnos.data(),       // material numbers
                                           names,               // material names
                                           1,                   // number of dimensions
                                           &nzones,             // pointer to dimensions
                                           0,                   // major order. row major is 0
                                           matlist,             // matlist
                                           mix_len,             // length of mix arrays
                                           mix_mat,             // mix_mat array
                                           mix_next,            // mix_next array
                                           NULL,                // mix_zone array (OPTIONAL)
                                           mix_vf,              // mix_vf array
                                           domain_name.c_str(), // domain name
                                           0);                  // allow mat0

        return mat;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetMaterial: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::GetSpecies
//
//  Purpose:
//      Gets the specset object from the Blueprint database and converts it 
//      into an avtSpecies.
//
//  Arguments:
//      domain     The domain of interest.
//      spec_name  The species set of interest.
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: Justin Privitera
//  Creation:   09/27/24
//
//  Modifications:
//
// ****************************************************************************
avtSpecies *
avtBlueprintFileFormat::GetSpecies(int domain,
                                   const char *spec_name)
{
    try
    {
        BP_PLUGIN_INFO("avtBlueprintFileFormat::GetSpecies " 
                        << domain << " "
                        << spec_name);

        Node n_specset;
        ReadBlueprintSpecset(domain,
                             spec_name,
                             n_specset);

        if (!n_specset.has_child("matset"))
        {
            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                 "specset " << spec_name << " is missing associated matset.");
        }

        const std::string matset_name = m_specset_info[spec_name]["full_material_name"].as_string();

        Node n_matset;
        ReadBlueprintMatset(domain,
                            matset_name,
                            n_matset);

        Node n_silo_matset;
        conduit::blueprint::mesh::matset::to_silo(n_matset,
                                                  n_silo_matset);

        Node n_silo_specset;
        conduit::blueprint::mesh::specset::to_silo(n_specset,
                                                   n_matset,
                                                   n_silo_specset);

        if (!n_silo_specset.has_child("speclist"))
        {
            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                 "intermediate silo representation of species set " << 
                                 spec_name << " is missing speclist.");
        }

        // first we need number of zones
        const int nzones = n_silo_specset["speclist"].dtype().number_of_elements();

        if (!m_specset_info.has_path(std::string(spec_name) + "/topo_name"))
        {
            BP_PLUGIN_EXCEPTION1(InvalidVariableException,
                                 "Species set information is missing topology for species set " << 
                                 spec_name);
        }

        Node n_tmp;
        n_silo_specset["nmatspec"].to_int_array(n_tmp["nmatspec"]);
        n_silo_specset["speclist"].to_int_array(n_tmp["speclist"]);
        n_silo_specset["mix_spec"].to_int_array(n_tmp["mix_spec"]);
        n_silo_specset["species_mf"].to_float_array(n_tmp["species_mf"]);

        const int    nmat        = n_silo_specset["nmat"].as_int();
        const int*   nmatspec    = n_tmp["nmatspec"].as_int_ptr();
        const int*   speclist    = n_tmp["speclist"].as_int_ptr();
        const int    mixlen      = n_silo_specset["mixlen"].as_int();
        const int*   mix_spec    = n_tmp["mix_spec"].as_int_ptr();
        const int    nspecies_mf = n_silo_specset["nspecies_mf"].value();
        const float* species_mf  = n_tmp["species_mf"].as_float_ptr();

        avtSpecies *spec = new avtSpecies(nmat,        // number of materials
                                          nmatspec,    // number of species associated with each material
                                          1,           // number of dimensions in the speclist array
                                          &nzones,     // array of length ndims that defines the shape of the speclist array
                                          speclist,    // indices into species_mf and mix_spec
                                          mixlen,      // length of mix_spec array
                                          mix_spec,    // array of length mixlen containing indices into the species_mf array
                                          nspecies_mf, // length of the species_mf array
                                          species_mf); // mass fractions of the matspecies in an array of length nspecies_mf
        
        return spec;
    }
    catch (conduit::Error &e)
    {
        std::ostringstream err_oss;
        err_oss <<  "Conduit Exception in Blueprint Plugin "
                    << "avtBlueprintFileFormat::GetSpecies: " << endl
                    << e.message();
        BP_PLUGIN_EXCEPTION1(VisItException, err_oss.str());
    }
}

// ****************************************************************************
//  Method: avtBlueprintFileFormat::RegisterDataSelections
//
//  Purpose:
//   Used to support avtResolutionSelection & capture the selected lod.
//
//  Arguments:
//     sels:    data selection list from the pipeline
//     applied: pipeline handshaking for handling data selections
//
//
// ****************************************************************************
void
avtBlueprintFileFormat::RegisterDataSelections(
                                const std::vector<avtDataSelection_p>& sels,
                                std::vector<bool>* applied)
{
    for(size_t i=0; i < sels.size(); ++i)
    {
        if(strcmp(sels[i]->GetType(), "avtResolutionSelection") == 0)
        {
            const avtResolutionSelection* sel =
                static_cast<const avtResolutionSelection*>(*sels[i]);
            this->m_selected_lod = sel->resolution();
            (*applied)[i] = true;
        }
    }
}
